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Abstract 

The first part of this thesis investigates the role of structured models in autonomous 
motor learning. Any autonomous system, such as the human motor system, has only 
the internal consistency of its various sensors to rely upon for model building (learn¬ 
ing). To study the possibility of learning structured models from internal consistency 
constraints, the specific problem of learning the kinematic parameters (relative link 
orientations and lengths) of general revolute joint manipulators is explored. First it 
is noted that a manipulator may form a mobile closed kinematic chain when inter¬ 
acting with the environment, if it is redundant with respect to the task degrees of 
freedom (DOFs) at the endpoint. Then it is demonstrated that if the mobile closed 
chain assumes a number of configurations, then loop consistency equations permit 
the manipulator and task kinematics to be calibrated simultaneously using only the 
joint angle readings; endpoint sensing is not required. Example tasks include a fixed 
endpoint (0 DOF task), the opening of a door (1 DOF task), and a point contact (3 
DOF task). Identifiability conditions are derived for these various tasks. The method 
is demonstrated for calibration of the Utah-MIT Dextrous Hand, and is generalized 
to hand-eye calibration. 

Part two focuses on the control of mechanical compliance during normal human 
elbow joint movement. In contrast to the first half of the thesis, this part stresses the 
experimental validation, rather than the formation, of control theories. Time-varying 
compliance estimates are made while subjects are executing normal movement. The 
estimates are made possible by the development of (1) a high performance wrist- 
mounted airjet thruster and (2) novel time-varying system identification techniques. 
The results indicate that the stiffness of the arm is low and is modulated during 
movement. The stiffness drops as soon as the movement starts and rises just before 
reaching a target. The implications of this and other findings are discussed in the 
context of feedforward control, compliance control, and equilibrium point control 
theories. Physiological mechanisms for stiffness modulation are also discussed. 
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Chapter 1 


Introduction 


The goal of this research is to understand the control of human arm movement. The 
view taken is that there are common problems involved in controlling any general 
purpose arm, whether it is biological or artificial. The intention is to uncover some 
of these problems, develop general methods for solving them, and infer the particular 
solutions taken by the human system. Ultimately, such solutions should reflect our 
motor competence, while adhering to the inherent mechanical and neural constraints. 

The investigation focuses on two aspects of the human movement control prob¬ 
lem: (1) learning kinematic sensorimotor transformations, and (2) controlling rapid 
arm movements. The first is investigated by implementing and evaluating calibration 
methods on robot manipulators. The goal is not so much to determine the details of 
human calibration, but instead, to gain an appreciation for the difficulties involved, 
and ultimately to derive appropriate experimental questions for future human stud¬ 
ies. The second topic is investigated by experimentally determining the elbow joint 
mechanical properties during rapid movement. These time-varying mechanical prop¬ 
erties are used to assess current theories of arm control. 


1.1 Background and Significance 

The control of arm movements is a deceptively difficult problem. Over the years 
researchers have come up with various ideas that are likely to be relevant to an 
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Figure 1-1: Feedback control of a joint angle trajectory 9. 

eventual theory of human motor control. A synopsis of some of these ideas provides 
a good starting point for the work contained in this thesis. The reader familiar with 
the context of this research may skip this section. It is intended to tie together the 
two loosely related parts of the thesis, and follows along the lines of [1] [2][3]. 

1.1.1 Feedback Control 

Given the complexity of both our motor system and the environment, it is reasonable 
to speculate that arm movements are made under strict feedback control (Figure 1- 
1). Take the simple task of extending the elbow joint. Muscle spindles, Golgi tendon 
organs, and joint receptors all contribute information about the actual trajectory of 
the elbow joint. Under linear feedback control the torque applied by muscle activation 
is made proportional (by a proportionality factor called the gain) to the difference 
between the actual and desired arm trajectory. Historically, Merton proposed the 
use of feedback control in controlling limb movements; he hypothesized that the 
gamma motor neurons are controlled and the main alpha motor neurons are driven 
from spindle feedback alone (see review [4]). Subsequent deafferentation studies and 
fusimotor neuron recordings proved him wrong, although some form of feedback is 
definitely present [4]. Before discussing the alternatives to feedback control, it is worth 
considering why feedback control alone cannot suffice to explain human performance. 
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Any feedback control scheme is by definition error driven. Thus, there are two 
ways to apply a greater force to speed up a movement: (1) tolerate greater errors, 
and (2) increase the feedback gains. 

Errors depend on the control variable 

The errors that are tolerated are task specific and depend upon the choice of 
the control variables. For example, in hammering in a nail the hand paths may 
be variable, while the impact location must be precisely controlled. In contrast, 
visually-guided pursuit requires continuous control of the hand trajectory (Sheridan 

[5] provides a good summary of this literature). Compromises between these two 
extremes have been observed for difficult pursuit-tracking tasks. For instance, Novas 

[6] found that when he tricked subjects by partially locking the target position to 
their hand movements (i.e., a “carrot on a stick” set-up) then movements were made 
intermittently. Bekey [7] also noticed evidence for discrete movements in pursuit 
tracking. He speculated that a visually-guided operator may be modeled as a sampled- 
data control system, sampling at 2-3 Hz (see also Crossman and Goodeve [8]). 
Stability, gain and phase 

As the loop gain of the feedback system exceeds unity instability may result. 
For example, consider the effect of the large neuronal propagation delays (spinal 
feedback delay exceeds 25 ms, and visual feedback delay exceeds 100 ms [9]). For 
a linear feedback control system pure delays can lead to instability as follows: if a 
frequency component of the input signal (e.g., desired position) has a period of twice 
the duration of the pure delay then the feedback will be 180 degrees out of phase - 
effectively providing positive feedback. Such oscillations will grow if the loop gain at 
that frequency is greater than one. 

Linear stability ideas generalize to non-linear systems such as the motor system. 
One approach is to generalize the notion of gain (using extended spaces [10]). Sta¬ 
bility is guaranteed if the product of the loop “gains” is less than unity (small gains 
theorem) [10]. Loop transformations make this theorem less conservative. A sec¬ 
ond approach is to generalize the linear systems result that keeping the open-loop 
phase shift between ±180 degrees ensures stability [11]. Stability is guaranteed for 
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the closed-loop system if both the forward path (e.g., the manipulator structure) and 
the reverse path (e.g., the feed-back controller, or the environment) are dissipative 
(strictly passive) [10]. This passivity theorem explains why simple proportional-plus- 
derivative feedback control of manipulators works at all (without gravity). 
Compliance is important 

High feedback gains also compromise the need for compliance. In the event that 
a manipulator runs into an unexpected object, the more compliant the manipulator 
is the less likely it is to break the object (or itself). The use of force feedback from 
the hand can improve the compliance of a high gain position controlled manipulator, 
but force sensor delays and non-colocated actuators lead to acute stability problems. 

In summary, using strict feedback control forces difficult trade-offs: speed, accu¬ 
racy and high gains versus stability, compliance (low gains), and energy dissipation 
(re passivity theorem). This is not to say that feedback control is not used by the 
motor system. On the contrary, reflexes are used, and in fact can be driven to in¬ 
stability (for example, with 8-12 Hz elbow joint perturbations; see Joyce and Rack 
[12][13]). The point is that the performance humans exhibit should be a result of 
more than feedback control alone. 

1.1.2 Pre-Planned Movements 

The problems with feedback control are countered by using pre-programmed (feedfor¬ 
ward) movements. See Figure 1-2. This was realized early on with the “motor tape” 
idea, and has taken on various theoretical forms in robotics [14] [15] [16] [17]. The cost 
of such an approach is the need for complex internal models. That is, the motor 
system and the environment must be represented, knowledge of the parameters of 
this representation must be learned, and appropriate feedforward trajectories must 
be computed. 

Of course, both feedback control and knowledge intensive planning can be com¬ 
bined. One method is to put knowledge of the non-linear arm dynamics into the 
feedback controller — as in model-referenced adaptive control [18]. Adaptive control 
has been used to model the details of the vestibular-ocular reflex (VOR) [19] (for a 
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Figure 1-2: Feedforward control of a joint angle trajectory 0. The feedforward control 
attempts to invert the arm dynamics V. 



Figure 1-3: Feedforward + feedback control of a joint angle trajectory 6. 

review of the various oculomotor control models see [20]). A feedback controller may 
also be used in cascade with a feedforward controller to reject unmodeled dynamics 
(e.g., [16]). See Figure 1-3. Owing to the non-linear inertial and gravitational link 
interactions - even at low speeds [21] - it is probable that some knowledge of the 
inertial properties of the arm are used for control. In human arm movements there 
is reason to believe that feedback is only used intermittently — to set-up successive 
ballistic movements, and avoid instability [8]. 

1.1.3 Mechanical Constraints Augment Active Control 

Another emerging theme in motor control research is an idea mechanical engineers 
sometimes use: it is advantageous to design a mechanical device such that its un¬ 
controlled passive behavior is close to the desired behavior. For example, consider 
McGeer’s [22] passive walking machine. Inspired by a child’s toy, McGeer noticed 
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that walking resembles the movement of a wheel, and thus designed a 2-legged device 
that emulates a wheel. Of course, actuators were needed to make the machine more 
versatile, but these actuators augmented, rather than ignored, the inherently stable 
passive walking behavior of the machine. In a similar manner, evolution of the the 
human motor system may have simplified the required control effort with mechanical 
design. 

In motor control and robotics numerous mechanical simplifications have been no¬ 
ticed: the human’s spherical wrist and shoulder joints allow potentially simple closed 
form inverse kinematic solutions [1]; the redundant arm geometry aids in singularity 
avoidance [23]; the use of contact friction helps in manipulation [24]; symmetry can 
simplify legged locomotion [25]; and finally, the mechanical properties of the muscles 
provide passive feedback control [8] [26] [27] [28]. 


1.2 Scope of Thesis 

The thesis concentrates on two specific problems, through which the above control 
issues are investigated. The first is the control of statically positioning the hand in 
space without visual guidance. This kinematic task is the most clear example of the 
use of feedforward control in the motor system. There must be some internal model 
(kinematic transformation) that maps hand locations into joint angles (or muscle 
lengths). A major issue and the focus of this thesis is the representation and learning 
of this kinematic transformation. 

The investigation of kinematic learning is organized as follows: Chapter 2 discusses 
the problem constraints and possible solutions. Chapter 3 presents a novel approach 
to learning the kinematics of the arm. Chapter 4 presents an implementation for 
calibration of the Utah-MIT Dextrous Hand. Chapter 5 generalizes Chapter 4 to 
calibration of a hand-eye system. Chapter 6 discusses the relevance of the methods 
presented in the previous chapters to human motor control. 

The second problem investigated is the control of rapid arm movement. This 
dynamic control problem has received considerable attention in both motor control 
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and robotics research; thus there axe ample theories of arm control. What is missing is 
experimental evidence to distinguish theories of feedforward versus feedback control, 
and to understand the role of intrinsic muscle properties. Thus, the second part of 
the thesis focuses on measuring the mechanical properties of the human elbow joint 
during movement. These properties are used to assess current theories of arm control. 

This study represents the beginning of a series of planned experiments that will 
eventually investigate whole arm movements. An important contribution is the de¬ 
velopment and verification of instrumentation and system identification methods re¬ 
quired for these future studies. 

Chapter 7 discusses the various movement control hypotheses and how mechani¬ 
cal impedance estimates may distinguish these hypotheses. Chapter 8 discusses the 
previous research with an emphasis on instrumentation methods. Chapter 9 discusses 
the methodology for the present studies. Chapter 10 presents the results from elbow 
joint movements. Chapter 11 assesses theories of motor control in light of Chapter 
10 . 
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Chapter 2 


Learning Sensorimotor 
Transformat ions 


The human motor control system is distinguished from current robotic control sys¬ 
tems by its remarkable ability to calibrate itself and adapt. Starting from clumsy 
infant movements, it is able to gradually develop itself to the point where walking, 
reaching, writing, playing tennis, etc., are all effortless. This adaptability entails 
autonomously acquiring knowledge (internal models) of the mechanical properties 
of the links, actuators, sensors, and environment. While current robotics systems 
may rely on an engineer to provide such knowledge, robotics systems are fast becom¬ 
ing sufficiently complex that human-like autonomous learning will soon be essential. 
Thus, both from the point of view of understanding the human motor control system, 
and of designing advanced robotic systems, it has become necessary to investigate 
autonomous learning. 

Any autonomous system, such as the human motor control system, has only the 
internal consistency of its various sensors to rely upon for model building (learning). 
To develop the notion of learning from internal consistency constraints, the problem 
of learning the kinematic parameters (relative link orientations and lengths) of gen¬ 
eral revolute joint manipulators is focused on. It is demonstrated that, provided the 
manipulator has sufficient redundancy, the kinematic parameters may be estimated 
using only the consistency constraints among the joint angle sensors. These consis- 
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tency constraints are generated by having the manipulator form a mobile closed-loop 
kinematic chain (e.g., as may be done by opposing one’s thumb to fore finger). The 
method has been implemented for calibrating the Utah-MIT Dextrous Hand, and as 
well extended to include autonomous camera calibration for hand-eye coordination. 


2.1 Terminology 

It is convenient to identify each joint sensor reading (e.g., joint angle, or muscle length) 
with a distance along a separate coordinate axis in a multi-dimensional space. Such 
an intrinsically defined space is referred to as a joint space. Likewise, the visual 
system has a intrinsic coordinate system that is referred to as retinotopic space. More 
generally, the term sensor space refers to other sensory modalities. We will also use 
the term to refer to the It is also useful to speak of a set of three-dimensional (3-D) 
world coordinates (e.g., Cartesian coordinates) that are fixed to the body, head or 
world. Such coordinates are often defined by the natural constraints of a task. This 
3-D coordinate system is variably referred to as task space, hand space, or world-based 
coordinates. See Figure 2-1. 

The transformations between the various sensory modalities considered in this 
part of the thesis are static, and are thus referred to as kinematic transformations 
- as opposed to dynamic transformations that convert muscle forces into joint angle 
trajectories. 

2.2 Problem Statement 

The relationship between our many sensors and muscles is non-linear and non-unique, 
yet we are able to utilize information gathered from one sensory modality (e.g., vision) 
to coordinate and control another (e.g., an arm) without continuous feedback. Our 
ability to make transformations between sensors is what underlies observed motor 
equivalence, such as our being able to write comparably with any appendage. As 
Held [29] points out, this ability to transform sensory information cannot be hard- 
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(a) (b) 

Figure 2-1: (a) Task space and (b) joint space. 

wired at birth: the transformations must adapt (calibrate) whenever growth, injury, 
or artificial manipulations occur. The central problem is to determine how calibration 
of these transformations is achieved. 


2.3 Representation 

The answer to the question of how calibration occurs depends on the representation 
of the transformations to be calibrated. The simplest representation is that of a large 
memory (look-up table) that relates the sensors in a discretized form (e.g., see [30]). 
For example, the transformation that maps joint angles of my left arm into joint 
angles of my right arm (when the hands are in the same location) can be represented 
by remembering all points in the right arm joint space that correspond to each point 
in the left arm joint space. This map is not one-to-one and requires substantial 
memory — perhaps even too much for our highly parallel brain architecture. On the 
related problem of determining the appropriate torques for moving the arm along a 
specified trajectory, there have been numerous proposals to reduce the memory size: 
hash tables [14], variable density table filling [17], and state-space methods [15]. 

On the other extreme, the transformation may be represented by structured mod¬ 
els with parameters that reflect the mechanical constraints of the sensorimotor system 
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(e.g., see [16]). This is the representation developed in this thesis. The hypothesis 
is that we have a priori knowledge of the structure of our limbs and sensors. For 
example, it is useful to assume: 

• rigid limbs 

• revolute joints 

• linear joint angle sensors 

These assumptions may be refined as needed (e.g., muscle geometry may be included), 
but such additional complexity would obscure, rather than clarify, the calibration 
problem and so is not introduced here. 

The above structural assumptions allow the development of parametric equations 
representing the sensory transformations. Calibration then proceeds by adjusting 
the parameters in the model (e.g., link lengths, and orientations) until the desired 
transform between any two sensory modalities is achieved. 

The benefits of this approach are that less memory is needed, and more generality 
and abstraction (see Section 3.4) is possible after calibration. This is because the 
structure of the model itself constrains the class of possible transformations a priori. 

A model-based approach also has a cost: (1) the kinematic parameter estimation 
usually involves a difficult non-linear search and (2) the structures of the models of our 
sensors and limbs are necessarily approximate; thus, accuracy is compromised. These 
two points are not separate; model accuracy may be exchanged for complexity of 
parameter estimation. For example, in camera calibration Direct Linear Transforms 
(DLT’s) [31] are used extensively for their ease of parameter estimation, yet the DLT 
method implicitly approximates the camera model with ideal projection properties, 
and it over-parameterizes the equations to make them linear. 

Finally, the use of structured models does not preclude memory-based representa¬ 
tions. Memory-based methods may be used locally to improve a structured model’s 
accuracy (trajectory learning) [16]; or even more extreme, they may be used to ob¬ 
viate the need to calculate in verse-kinematics. These issues are tangent to the main 
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topic of calibration. In calibrating sensor transformations there is still a more subtle 
problem that must now be addressed. 


2.4 Learning Transformations to a 3-D World 

Learning transformations to task space is more difficult than learning transformations 
between sensors. The difficulty is that there are not sensors that directly provide 3-D 
coordinate information (or some metrical equivalent). It is not sufficient to assume 
that the visual system provides task-space coordinates, for the visual system must 
itself be calibrated. Thus, the methods of learning by association of input-output 
pairs must be re-thought - regardless of how the transformations are represented. 

2.4.1 Is there a Task Space? 

The problem of learning to transform sensor data into a common abstract reference 
frame (e.g., task space) is so serious that it is tempting to think that we might not 
maintain such a 3-D world representation. As discussed in Chapter 6 this thought 
is probably wrong. It would be hard to imagine doing the type of reasoning that 
we are capable of without an interned 3-D representation of the world; without it, we 
would have to reason in joint space (or retinotopic space) where distances in the world 
are not preserved, rigid objects distort as they are sensed in different locations, and 
kinematic redundancy turns points in the world into surfaces. Thus, short of claiming 
that we do not maintain a metrical representation of the world (instead using only 
topological properties that are preserved by the sensory transformations [32, page 
49]), it must be concluded that we somehow learn transformations from our sensor 
spaces to a 3-D task space. 

2.4.2 Teacher Based Methods Do Not Suffice 

Assuming that we do learn transformations to a 3-D task space, and assuming that 
task space is distinct from joint space or retinotopic space, it is possible to re-assess 
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the old dispute between memory-based and structured model-based learning - at least 
for kinematics. 

To learn the transformations from task space to joint space some external agent 
(teacher) must provide the task-space coordinates of the hand location corresponding 
to a given arm configuration (point in joint space). This requires the equivalent of 
the construction of a measurement machine against which to calibrate our sensors. 
Although this is a standard approach in robotics, it is an unreasonable proposition 
to explain how an infant learns to transform sensory data into a common task-space 
coordinate frame. Thus, black-box learning methods, such as unstructured neural 
networks or memory look-up tables, do not suffice for learning to transform data to 
an abstract task space. 

2.4.3 A Priori Knowledge of Physical Structure Enables 
Abstraction 

The situation is more promising for structured models. As mentioned, structured 
models have a representation of the 3-D world built into them. They provide the 
necessary constraints by which to compare disparate sensory information. Calibra¬ 
tion may proceed by adjusting the parameters of the transformation models so as to 
maximize the consistency among the sensors. Once the parameters are estimated, 
transformations to an abstract task space can be computed analytically. Let us see 
how this might work in a very simple example. 


2.5 Simple Planar Calibration Example 

Consider a 3-D OF planar manipulator making a point contact with the ground. This 
manipulator is redundant with respect to the point contact constraint, and thus forms 
a mobile 4-bar closed linkage. See Figure 2-2. (Figure 2-2 may be viewed alternately 
as a 1-DOF manipulator being tracked by a 2-DOF manipulator.) The ground is 
considered to be the fourth link with a fixed length of a 4 . The goal is to determine 
the kinematic parameters aj, ct 2 > <* 3 * and a 4 from the joint angle readings $i, $2 and 63. 
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Figure 2-2: A 3-DOF planar manipulator with a fixed point contact to ground. 

The three kinematic constraint equations (at the ith configuration) may be written 
as: 


axcos^i) + a 2 cos(0\ + 0' 2 ) + o 3 cos(^ + + 0 3 ) + a 4 = 0 ( 2 - 1 ) 

0 !sin(^j) + a 2 sin{0\ + $* 2 ) + a 3 sin($\ + + ^ 3 ) = 0 ( 2 - 2 ) 

E*j = 0 (2.3) 

j=i 

Evidently, the length parameters may be scaled arbitrarily and still satisfy equations 
(2.1) and ( 2.2). For this reason one link length must be defined as unity. Let a 4 = 1 . 
To continue, each additional configuration of the mechanism provides two additional 
position equations. Placing the equations from two configurations into matrix form 
we have: 


Q. = C<£+y 


(2.4) 
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where <p = (a x a 2 a 3 ) T , y = ( 1 0 1 0) T , ^*- — 9\ and 


cos(<f>\) cos(<f> 2 ) cos(<f> 3 ) 
sin($) sin(^) sin(<f>\) 
cos(<f>\) cos{(j>\ ) cos(^> 3 ) 
sin{(j>\) sin(<f>\) sin(<f >\) 


(2.5) 


The least squares solution is £ = (C T C)~ 1 C T (—y). A unique solution is guaranteed 
provided that the columns of C are independent. Observe that the columns of C 
will only be dependent if the mechanism happens to be a parallelogram (that is, 
4>\ = $,). Once calibration is complete (i.e., <£_ is determined), arbitrary positioning 
of the manipulator in task space is possible with the kinematic equations ( 2 . 1 ) and 
( 2 . 2 ). 


2.6 Closed-Loop Kinematic Calibration 

In the next three chapters the above example will be generalized, and referred to 
as closed-loop kinematic calibration. In contrast to the planar example, the general 
calibration equations are non-linear. Thus, the calibration must proceed iteratively, 
starting with initial parameter estimates, and may arrive at local or multiple solutions. 
It is crucial that eventually the parameters are uniquely determined by the data and 
constraint equations; if they are not, they may only model the training data set 
locally, and may not be of use in computing transformations to arbitrary points in 
task space. 

Though the calibration method will be developed strictly in the context of robotics, 
the generalization to motor control should be clear from the discussion in this chapter. 
After developing the method, Chapter 6 will further discuss the relevance of this 
model-based approach to motor control. 
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Chapter 3 


Closed-Loop Kinematic 
Calibration 


3.1 Introduction 

Kinematic calibration is important for model-based robot control [16]. We [33] and 
many other authors (see review paper [34]) have developed open-loop methods that 
estimate the geometric and static non-geometric kinematic parameters of open-chain 
manipulators, by relying on special purpose pre-calibrated endpoint locating systems, 
such as precision points or camera-based measurement systems. Section 3.2 summa¬ 
rizes our open-loop method, which is the starting point for our new method; new 
results on the identifiability of the open-loop method are also provided. 

Our new method, which we call closed-loop kinematic calibration, eliminates the 
need for endpoint locating systems: if a manipulator is formed into a mobile closed 
kinematic chain, then its joint angle readings alone are enough to identify the kine¬ 
matic parameters. A manipulator may form a mobile closed-loop kinematic chain if it 
is redundant with respect to its endpoint constraint (task). Section 3.3 considers the 
simplest endpoint constraint: the position and orientation of the endpoint are fixed 
relative to the base link. For this O-DOF task, the manipulator must be redundant 
(> 7 DOFs) to form a mobile closed loop. For each loop configuration there are three 
position and three orientation loop consistency equations. If the closed loop is placed 


28 



into n configurations (with the same endpoint location), there result 6 n equations 
that may be solved for the unknown parameters. 

An equivalent scenario is two manipulators rigidly attached together at their end¬ 
points with a combined total of DOFs > 7. The last link of the second manipulator 
may be defined as the base, and the entire closed kinematic chain may be viewed as 
a single equivalent manipulator. 

Section 3.4 then considers two endpoint constraints with passive degrees of free¬ 
dom: 1) rotation about a passive revolute joint (e.g., opening a door), and 2) rotation 
about a passive spherical joint (e.g., point contact). For constraint 1, a non-redundant 
6-DOF manipulator can form a mobile closed loop, while for constraint 2 a 4-DOF 
manipulator could suffice. In general, for each passive degree of freedom introduced 
into the endpoint constraint one less degree of freedom is required by the manipula¬ 
tor. As a corollary, the geometry of the task can also be identified, for example, the 
position and orientation of the revolute joint in constraint 1. 

Several technicalities were overcome in developing this closed-loop method. A the¬ 
orem was developed to determine which parameters are identifiable in the consistency 
equations. We also show how the passive DOFs can be eliminated from the endpoint 
constraint for the two cases studied, and mention how to do it in general. Thirdly, we 
apply a Newton-like search method for the kinematic parameters, which starts from 
an initial guess at the parameters. Simulations will demonstrate the convergence of 
the method. Finally, the manipulator must be able to make constrained internal joint 
movements, without knowing the true kinematic parameters or producing excessive 
internal or endpoint forces. 

Closed-loop kinematic calibration is related to mechanism synthesis [35], which 
characterizes closed-loop mechanisms with one degree of mobility through relative dis¬ 
placements of designated input and output angles. By eliminating (with difficulty) 
the five unspecified DOFs from the kinematic equations, a displacement equation 
results that is a 16th order polynomial in the tangent half-angles of the input and 
output angles [36]. A difference from mechanism synthesis is that serial-chain ma¬ 
nipulators typically have sensors on all the joints, and so eliminating the unsensed 
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Figure 3-1: A single closed-loop kinematic chain formed by a redundant manipulator 
or by dual manipulators. 


passive DOFs at the endpoint from the kinematic equations is considerably easier. 
Portions of this work have been previously reported [37], [38], [39]. 


3.2 Open-loop kinematic calibration 

This section presents our method for open-loop kinematic calibration [33], which 
serves to set the basic concepts and mathematics from which the closed-loop method 
is derived. New results in identifiability are presented for our open-loop method; these 
results apply more generally to similar methods that have appeared in the literature 
[34]. 


3.2.1 The Kinematic Model 

Both geometric and non-geometric parameters are required for kinematic calibration. 
The Denavit-Hartenberg (D-H) convention [40] is employed for the geometric param¬ 
eters (Figure 3-2). For a manipulator with n DOFs, the end effector is located by the 
position vector p’ and the orientation matrix R’: 


n 

^ 1 + fljXj 

3 =1 

n 

(3.1) 

n R.(«j)R.(a‘) 

j=l 

(3.2) 
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Figure 3-2: Denavit-Hartenberg coordinates and tip vector b^. 

where R 2 (</>) and R*(^) are 3x3 rotation matrices about z and x axes by the angle 
<f >, and the subscript c indicates that the position and orientation are computed from 
the model. The superscript i refers to the configuration of the manipulator, since in 
kinematic calibration it is placed into a number of configurations 0’ = (0 \,..., 0^), 
i — 1 The required geometric parameters are Sj, ocj and aj, for links j = 

1 ,... ,n. 

The non-geometric parameters are focused at a joint, and reflect errors between 
the true and measured joint angle; sources of error include backlash, gear eccentricity, 
joint compliance, and joint angle offset. We model only the joint angle offset error 
6°/*, which needs to be identified. It is related to the actual 0'- and measured Oj D-H 
joint angles by 0'- = Oj + 0°j*. All of the unknown kinematic parameters are placed 
into a single vector <£_ = (0 o// , 3, a, a), where s = (a 1? ... ,s n ), etc. 

Instead of the orientation matrix R*, it is convenient to represent the orientation 
by the vector r’ = (<^ x , <f> y , representing the xyz Euler angles. The computed 
endpoint location a£ = (p*, r*) may then be written as: 

a£ = Z(ff,£) (3.3) 

where the function / is derived from (3.1)-(3.2). It’s exact form is not required here. 
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3.2.2 Iterative Identification 


To estimate p, the manipulator must be moved into an adequate number m of con¬ 
figurations, in consideration of the large number of parameters in p and of statistical 
averaging. At each configuration i the actual endpoint location a£ is measured. The 
goal is to determine the p that best predict from the kinematic model (3.3) all of the 
endpoint measurements X = (a£,..., x£*): 

X = Ht) (3.4) 

where ?■(£) = (/(«’,£),...,/(« m ,£)). 

Solving for p_ from (3.4) is a nonlinear estimation problem, which can be done by 
linearization and iteration: 


AX = CAp (3.5) 

where C = dT/dp. We can consider AX = (A* 1 ,..., Ax m ), with Ax* — a£ — a£, as 
the location errors. Similarly, Ap = p — p f) is the error in the total parameter set, 
where p^ is the current estimate and p is the corrected estimate. In Ap, As = s — s^, 
etc. An estimate of the parameter errors is provided by minimizing LS = (AX — 
Cp) T (AX — Cp), which yields 


Ap = (C r C)~ 1 C r AA' (3.6) 

Finally, the guess at the parameters is updated as p = p^ + A p and the iteration 
continues until AX —► 0. 

The basis for linearization is the assumption that a£ is close to a£. Then 

A£ = - a£ = (dx% dy\ dz\ 8x\ dy\ dz { ) (3.7) 

where Ap* = (dx *, dy *, dz*) is the incremental position error and Ar* = (dx*,dy*,dz*) 
is the incremental orientation error. When p^ is far from the final values, problems 
with this approach may occur, as discussed later. 
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The Jacobian C is often formed through the use of differential homogeneous trans¬ 
formations, but to address identifiability we find the use of screw coordinates more 
convenient and transparent [41], [42], [43]. The differential of (3.3) is 

A *' = -g£ Ai + H** + ^ + W^ (3 ' 8) 

where 

Qj ^ 

5^ 5£ 5s 5a da 

and 


C 


c 1 

c n 


(3.9) 


Subsequently, we use the abbreviation = dx*./d/3, as each matrix in (3.8) is a 
Jacobian with respect to a particular parameter; thus 3g is the ordinary Jacobian 
related to joint angle displacement. 

The endpoint variation A®’ can be considered an instantaneous screw displace¬ 
ment composed of an incremental translation Ap’ and rotation Ar‘, caused by the 
combined variation in all of the parameters. For example, a variation A Sj contributes 
As_,z’_i to Ap‘. A variation A aj contributes (Aa_,)xj to Ar’ and (A<*j)x) x b* +1 to 
Ap’, where b* +1 is a vector from the jth coordinate system to the endpoint (Figure 
3-2). The Oj and aj parameters are treated analogously. In total, 


Ap*' = *5-i x b 5 A % + z 5-i As ; + x 5 x b 5+i Aa i+ x 5 Aa i ( 3 - 10 ) 

j =i 

n 

Ari = E z 5-v+ x 5 Aa i (3-n) 

3=1 

Comparing to (3.8), it is seen that 
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(Ji)i = 

X*- 

j 

» = 

z 5-i 

II 

*5-i x b i 

, rt)i = 

x 5 x b 5+i 


0 


0 


. z 5-i . 


11 i 


(3.12) 

where the yth column of each Jacobian is indicated by the subscript, and represents 
a screw coordinate for variation in the parameter f3j. 


3.2.3 Identifiability 

Next we derive several results pertaining to the identifiability of <p in (3.4). First we 
establish that the solution cannot be globally unique. 

Theorem 1 There are at least 2 n_1 solutions <£_ to (3.4). 

Proof. We presume that at least one solution <£_ exists, because the data come from 
a physical system. Additional solutions may be derived from this <p. There are two 
possible parameter sets per joint. For a fixed zj, the x*- axis can be made to point in 
opposite directions by adding 180° onto to accommodate this change, the sign is 
changed on aj and aj while Sj is unchanged. At the endpoint, the directions xj, and 
z' n are specified by the position requirement. Hence we have generated 2 n_1 solutions 
from the original solution <£. □ 

Though there are multiple solutions, in practice kinematic calibration starts off 
with a rough estimate and searches locally for a solution. Thus, the relevant 
question is whether or not there is a unique solution within a small region of the 
parameter space. We draw on some results from differential topology [44], [45]. 

Definition 1 (Locally unique) A solution <£) to (3.4) is locally unique if it is the 
only solution in an arbitrarily small neighborhood (ball) around <p'. 

In addition to uniqueness, we also want to know if a well-behaved inverse function 
J-~ x exists to generate the solution pf. 
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Definition 2 (Locally identifiable) A smooth function X = iF(<p) is locally iden¬ 
tifiable at tpf if (1) is locally unique, and (2) there exists a smooth inverse function 
J-~ x such that .F -1 (.F(^)) = <p. 

A smooth function has continuous partial derivatives of all orders. Notice that the 
existence of a smooth inverse function in part 2 of the definition does not guarantee 
uniqueness (part 1). 

If the number of equations equals the number of unknown parameters, then local 
identifiability is equivalent to requiring that T is a local diffeomorphism. It is es¬ 
tablished in [44] that J- is a local diffeomorphism if and only if the square Jacobian 
matrix C = \dT / d<pf\ is non-singular. This motivates the following results for when C 
is not necessarily square. Lemma 1 is a general result that applies to all equations of 
the form (3.4). Lemmas 2 and 3 are particular to the kinematic calibration problem. 

Lemma 1 Let be a solution to (3.4)- The Jacobian C has full rank if and only if 
the parameters <£^ are locally identifiable. 

Proof. Assume C has full rank. Let <p be another solution to (3.4). A Taylor series 
expansion relates <£_ to <£f: 


T{<£) = T{<£) + C{<£ -£') + ••• (3.13) 

where F(<f) = .T’X<£>'). If <p is in a sufficiently small neighborhood of then the higher 
order terms after the first differential may be neglected. Then C(<p — ip') = 0. Since C 
does not have a null space, then <£_ = (£?, and <£f is locally unique. Furthermore, there 
exists a smooth inverse function J : 1 , as (3.6) suffices. 

Conversely, assume the parameters <£' are locally identifiable. Let T~ x be a smooth 
inverse function. Differentiating T~ 1 (T(<£)) = <£_ by <p yields [dT~ x /dX][dT/dp] = I. 
If C does not have full rank, then there exists a vector element n of its null space. 
Postmultiplying by n yields 0 = n, a contradiction. Hence C has full rank. □ 

If Jacobian C has dependent columns, the solution to (3.4) is not necessarily locally 
unique because the higher order terms of the Taylor series expansion of J~ may be 
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non-zero when evaluated at an element of the null space of C. The following two 
lemmas require the definition of C by (3.8); they do not hold in Section 3.3. 

Lemma 2 The Jacobian C does not have full rank if and only if there exists a constant 
linear relation among the x*- and z*- axes, that is, 

n 

0 = 5>; z 5-i + ^ x i (3.14) 

3=1 

for some constants Cj and kj, not all zero, for all configurations i — 1,..., m. 

Proof. Part 1: Assume C does not have full rank. Then there is a constant linear 
relation among the screw coordinates (3.12). Two cases must be considered. First, 
if this linear relation only includes the left two screw coordinates in (3.12), then the 
same linear relation must hold in the positional component of the screws, the first 
three rows that contain just the x*- and z‘- axes. Thus (3.14) holds for some constants 
Cj and kj, not all zero, for all configurations i = l,...,m. Second, if this linear 
relation includes the right screws in (3.12), then this linear relation must hold in the 
rotational component of the screws, the last three rows that also contain just the x) 
and z’ axes. Again, (3.14) holds. 

Part 2: Conversely, assume (3.14) holds. Then at least the screw coordinates (J a )j 
and (J 4 )j have a constant linear dependence, and C does not have full rank. □ 

Lemma 3 C has full rank if and only if cp' is locally unique. 

Proof. Assume C has full rank. By Lemma 1, t£? is locally identifiable and hence is 
locally unique. Next assume <£ is locally unique. If C did not have full rank, then 
by Lemma 2 we could add (3.14) to (3.1) to change the length parameters without 
affecting p*. The kinematic length parameters would not be unique, contradicting 
the assumption. Hence C has full rank. □ 

Theorem 2 (Identifiability) The parameters <£ are locally unidentifiable or locally 
non-unique if and only if there exist constants Cj and kj, not all zero, such that 
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(3.15) 


0 = E c d-t + MJ 

3 =1 

/or all configurations i = 1,... ,m. 

Proo/. This theorem follows from Lemmas 1, 2 and 3. □ 

Various categories of singularities (occurrences of (3.15)), which are generic to the 
closed-loop case as well (see Section 3.3), will now be enumerated. These categories 
are meant to be illustrative rather than exhaustive. Although not expressly discussed 
in each category, the real problem is associated with near singular situations, which 
cause intractable numerical sensitivity problems while solving for the parameters. 

Singularity 1: coordinate description. In the D-H convention, when there are 
two consecutive parallel joint axes, there is no unique common normal (Figure 3-2). 
Parallel axes imply: 

z 5 - z i-i = 0 (3.16) 

Thus (3.15) is true and the corresponding s,- and s, +1 may not be identified alone 
(although the difference Si — s l+1 can be identified). This situation can be avoided 
by changing the coordinate description of the parallel axes to a convention such as 
Hayati’s [46], which however may not be used exclusively because it too suffers from 
a parameter ambiguity when two consecutive joint axes are perpendicular. 

A revolute joint axis is a line vector, which is located by 4 parameters. Hence 
any coordinate description with greater than 4 parameters per link is singular (unless 
extra constraints are imposed). Similarly, a prismatic (linear) joint axis is a free vector 
defined by only 2 orientation parameters, and coordinate descriptions with more than 
2 parameters are singular. 

Singularity 2: insufficient excitation. If the mechanism is not moved into a suf¬ 
ficient number of configurations, then the data are not sufficiently exciting [47]. A 
small variation in each parameter of <£_ should cause an observable displacement in 
the end effector. The optimal data set would maximize the observable model error 
over variations in all of the parameters [48]. An impoverished data set would not 
be able to distinguish changes in particular parameters, which could vary arbitrarily. 
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A trivial example is an immobile joint whose axes x*-, z’-, x*_ x , and z*_ x are always 
linearly dependent. 

Singularity 3: transient singularities. During the course of the iterative search, an 
intermediate singular parameter set may be found even though the real mechanism 
may not have singularities. Simulations show that this situation is surprisingly com¬ 
mon when the initial guess is not close to the true solution. Since this singularity 
is associated with the algorithm, it may be avoided by the modified minimization 
criteria LS' = LS + AA^ r A <p. In addition to minimizing the end effector tracking 
error, LS' minimizes the variation in parameters so that at a potential singularity the 
arbitrary parameters tend to remain fixed. Minimizing LS' yields 

A<£ = (C T C + Xiy^AX (3.17) 

Iteratively applying (3.17) results in the Levenberg-Marquardt algorithm [49], [50]. 
The free parameter A determines the trade-off between a straight Newton iteration 
and a much slower gradient descent. 


3.3 Closed-Loop Kinematic Calibration 

We consider next a redundant manipulator (> 7 DOFs) rigidly attached to the ground 
at its endpoint. In general, the resulting closed-loop kinematic chain is mobile, since 
the fixed endpoint constrains only 6 of the DOFs of the manipulator. Our closed-loop 
method makes use of this mobility to kinematically calibrate the manipulator without 
endpoint sensing. The following observation is the key: the origin of coordinates can 
be placed at the fixed endpoint location and can be defined to have zero orientation 
and position. Hence a£ = 0, and no measurements are required because the actual 
endpoint location is known and is zero (by definition). 

Figure 3-3 illustrates this origin placement for a 7-DOF manipulator. The com¬ 
bined manipulator tool and ground may be viewed as a single effective link that 
connects the last to first joint of the manipulator. The end effector axes z\ and x^ are 
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made coincident with the base coordinate axes zj, and x l 0 . Note that the kinematic 
parameters of the ground link must now be identified. 

The mathematical development in the previous section then applies, with one 
modification. As before, the mobile closed chain is moved into a number of configu¬ 
rations. At each configuration i the endpoint error follows simply from (3.7): 

A®’ = = ~{dx \, dy ' c , dz' c , dx' c , dy[, dz\) (3.18) 

where Ap’ = (dx\, dy x c , dz' c ) is the computed position and Ar* = (dx]., dy' c , dz' c ) is the 
computed orientation. The iterative estimation procedure cannot be applied further 
without modification because the Jacobian C is singular. The position equations for 
the closed loop are: 

n 

Pc = £ = 0 (3.19) 

3 =1 

Hence the length parameters are linearly dependent and (3.14) is satisfied. Intuitively, 
the size of the manipulator can be scaled arbitrarily and still satisfy the loop closure 
equations. 

To proceed with our closed-loop method, it is necessary to specify one length 
parameter to scale the size of the mechanism. For example, suppose we set a\ = — 1 . 
We redefine a = ( 0-2 >•••,«») and remove the first column from the Jacobian J*, 
which redefines C 1 in (3.8). We may then proceed as before with the parameter 
identification. Similarly, if another parameter such as S 3 had been specified instead, 
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analogous changes in the definition of s and C* would be required. In the remainder 
of this paper, we vary which parameter sets the scale for convenience. 

3.3.1 Identifiability 

For the purposes of identifiability, we will proceed with the scaling = —1. The 
following identifiability results are couched in terms of the choice of a*, but a different 
choice would result in trivial changes to the results. Of course, it is necessary in the 
actual mechanism that a\ 0. We do not consider mechanisms that have all length 
parameters zero (i.e., spherical mechanisms). 

Theorem 1 and Lemma 1 apply intact to the closed-loop case, but Lemmas 2-3 
an d Theorem 2 require a slight modification. Redefine the endpoint position as 

r» n 

Pc = *i = E s ri-i + E (3.20) 

3 =1 3=2 

Also, we must make the following definition to treat a possible exception arising 
from x\ being in the expression for (J^) x . 

Definition 3 (Type-E Special Mechanism) A single loop closed kinematic chain 
is type-e if its screw coordinates in (3.12) satisfy the following exceptional relation: 

0 = E'ifJiJj+OjfJDi + ftWj + otJiJj (3.21) 

3=1 

for all configurations i = 1,..., m. qj and rj are arbitrary constants. 

Perhaps (3.21) never occurs, but for completeness we include its possibility. 

The new identifiability theorem can now be stated. 

Theorem 3 (Fixed Endpoint Identifiability) For a redundant manipulator with 
fixed endpoint forming a closed kinematic chain that is not type-e, the parameters 
are locally unidentifiable, or locally non-unique if and only if there exist constants Cj 
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(3.22) 


and kj, not all zero and not all Cj = Sj and kj = aj , such that 

o = E c ; z i-i + 

j=i 

for all configurations i = 1 ,... ,m. Furthermore, the parameters of a type-e mecha¬ 
nism are non-unique, if (3.22) holds for the Cj and kj restricted as above. 

Proof. With the added restriction that not all Cj = Sj and kj = aj Lemma 2 holds if Xj 
is eliminated from (3.14) by substitution of (3.20), provided that the mechanism is not 
type-e. Lemma 3 holds if its proof is modified to use (3.14), with x^ eliminated, and 
(3.20) instead of (3.1). Alternatively, for all mechanism types local non-uniqueness 
follows directly from (3.22) by adding (3.22) to (3.20) and remarking, as in the proof 
to Lemma 3, that the length parameters are non-unique. This theorem then follows. 
□ 

For type-e mechanisms the angle-dependent screw coordinates (J„)^ and (J^)j in 
(3.12) are linearly related by the length parameters as defined above. The orientation 
component of this linear screw relation (the last 3 rows) gives (3.14) with Cj = Sj and 
kj = aj for all joints j — 1,... ,n. This provides no addition information, as it is the 
same as (3.20), and case 2 in part 1 of the proof to Lemma 2 cannot proceed. We do 
not know if a type-e mechanism can actually occur and do not consider it further. 

We now discuss two additional singularities in the closed-loop calibration proce¬ 
dure to the three singularities in the open-loop procedure that also apply here. 

Singularity 4- Inherent singularities in the mechanism. Certain mechanisms have 
particular symmetries that allow the kinematics to be described in less than four 
parameters per joint. It is difficult to provide a general rule for when this will happen, 
but it is usually restricted to mechanisms of mobility one. A simple example is a 3- 
DOF planar manipulator that makes a point contact with the ground (Figure 2-2). 
If the resulting closed-loop, four-bar linkage happens to be a parallelogram, then the 
opposite x axes are always parallel: 


x 2 + x 4 = 0 


(3.23) 
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This satisfies the condition (3.22), and thus the lengths of the opposite sides, a 2 and 
a 4 , axe not identifiable (except as a sum). Clearly, this problem may be eliminated by 
having the manipulator change its endpoint location so that a parallelogram is not 
formed. 

Singularity 5: structural immobility . If a particular joint j is immobile, then 
two consecutive joint coordinates are fixed relative to one another. This implies 
that xj, Zj, xj_i, and Zj_ x have a constant linear relation (satisfying (3.22)), as 
these four vectors span a three dimensional space. Of course, it is not surprising 
that the parameters of the links connected by the immobile joint are unidentifiable. 
Conceivably a fictitious link that combined the two links could be defined and the 
rest of the mechanism identified. 

To proceed, it is necessary to spot immobile joints. Following the approach in 
[42], we first determine whether the mechanism is totally immobile. Since the classical 
mobility definition [35] does not suffice for special mechanisms, the following condition 
is derived. 

Lemma 4 A single-loop closed kinematic chain is mobile if and only if the columns 
of the Jacobian JJ are linearly dependent. 

Proof Let the screw coordinate $’• represent the jth column of J^. Since J*!?’ = fi, 
then 


0 = + • • • + (3.24) 

All Oj’s will be identically zero if and only if the joint screws $’• are independent. □ 
Next we determine whether a single joint, say joint 1, is immobile. From (3.24), 
link one’s instantaneous movement is = — 0^$^ — ••• — 0\% x 2 . For link one to 
move, Sj must be a linear combination of the other screws. Stated otherwise, the 
space span[$j] must intersect the space span^,..., $ 2 ] spanned by the other screws. 
The following result from linear algebra is useful [51]. 
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Figure 3-4: Anthropomorphic arm screw coordinate assignment. 


Joint 

& ( m ) 

Sk (m) 

SL (rad) 

i 0,t (rad) 

1 

1.694 

0.837 

3.774 


2 

1.622 

-0.627 

-0.553 


3 

1.000 

-0.100 

0.930 

■jTTjTM 

4 

-0.430 

-0.430 

2.040 


5 

0.540 

-0.600 

-0.150 


6 

-0.560 

-0.550 

1.350 


7 

-1.693 

0.711 

-1.859 



Table 3.1: 7 DOF mechanism: initial parameters. 

Lemma 5 Let A and B be subspaces of a vector space V such that V = A-\-B, where 
A + B = [v \ v = a + b,a £ A,b (E B\. Then 

K(A r\B) = K(A) + K(B) - K(V) (3.25) 

where K(W ) denotes the dimension of a vector space W . 

For the following lemma, identify A with span[$j] and B with span[$^,..., $j]. 
Lemma 0 Joint one will be mobile if and only if 

1 + Kispan ^,$*]) - Kispan^i ,..., $‘]) > 0 (3.26) 

Any joint’s mobility may be ascertained by (3.26) with the appropriate re-numbering 
of the links. 

As an example, consider the mechanism formed by rigidly fixing the hand of an 
anthropomorphic arm [23] relative to its shoulder (that is, imagine holding onto the 
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Joint 

s (m) 

a (m) 

a (rad) 

(rad) 

1 

1.594 

0.737 

3.604 

0.000 

2 

1.722 

-0.527 

-0.503 

0.000 

3 

1.000 

0.000 

0.530 

0.000 

4 

-0.330 

-0.330 

2.300 

0.000 

5 

0.440 

-0.440 

-0.900 

0.000 

6 

-0.660 

-0.550 

1.200 

0.000 

7 

-1.793 

0.911 

-1.459 

1.825 


Table 3.2: 7 DOF mechanism: calibrated parameters. 


Joint 

s (m) 

a (m) 

a (rad) 

(rad) 

1 

-1.694 

-0.837 

3.504 

0.100 

2 

-1.822 

0.627 

-0.553 

0.050 

3 

1.000 

0.100 

0.580 

0.070 

4 

0.430 

0.430 

2.380 

1.070 

5 

-0.540 

0.540 

-0.980 

0.080 

6 

0.760 

0.650 

1.280 

0.050 

7 

1.200 

0.760 

-1.390 

0.100 

8 

1.200 

1.500 

3.800 

0.100 

9 

0.600 

-1.400 

1.550 

0.200 

10 

0.300 

1.100 

-1.380 

0.300 

11 

1.300 

0.600 

0.880 

0.900 

12 

-1.982 

-1.839 

1.772 

0.400 


Table 3.3: 12 DOF mechanism: initial parameters. 

desk in front of you). Although this mechanism has a classical mobility of 1 (7 minus 
6), the elbow joint can be shown to be immobile. Consider the upper arm as the 
base link, so that $\ is the screw coordinate for the elbow (Figure 3-4). Without 
loss of generality, one of the three wrist joint axes may be defined to intersect the 
shoulder joint. Thus the three shoulder joint axes and this wrist joint axis are linearly 
dependent, and K(span [$' 2 ,..., $^]) = 5 whereas K (span [$\,..., $^]) = 6. Therefore 
(3.26) shows that the elbow joint is immobile. The solution to this problem is to 
relax the endpoint constraint so that the elbow is mobile; for example, only maintain 
a point contact with the ground, allowing arbitrary orientation of the hand. This 
solution requires a reformulation of the identification equations and is taken up in 
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Joint 

s (m) 

a (m) 

a (rad) 

9°” (rad) 

1 

-1.594 

-0.737 

3.604 

0.000 

2 

-1.722 

0.527 

-0.503 

0.000 

3 

1.000 

0.000 

0.530 

0.000 

4 

0.330 

0.330 

2.300 

0.000 

5 

-0.440 

0.440 

-0.900 

0.000 

6 

0.660 

0.550 

1.200 

0.000 

7 

1.100 

0.660 

-1.300 

0.000 

8 

1.100 

1.400 

3.900 

0.000 

9 

0.500 

-1.300 

1.400 

0.000 

10 

0.200 

1.000 

-1.300 

0.000 

11 

1.200 

0.500 

0.800 

0.000 

12 

-1.882 

-1.939 

1.722 

0.000 


Table 3.4: 12 DOF mechanism: calibrated parameters. 


Section 3.4. 

3.3.2 Simulations 

This section presents two simulations, one for a 7-DOF manipulator and the other 
for two 6-DOF manipulators rigidly attached at their endpoints. In these and all 
subsequent examples the rank of the matrix C was monitored to avoid singularities. 
Further, C has full rank for the actual parameters, and thus all mechanisms are 
identifiable. 

Example 1. A 7-DOF manipulator, with actual D-H parameters in Table 3.2, was 
formed into a closed loop mechanism by the end effector grasping the ground at a fixed 
arbitrary position. This mechanism was simulated in 40 distinct configurations = 
0 to 0.5 rad). Starting with the initial guesses in Table 3.1 and with the definition 
S 3 = 1, the simulated joint angles were fed into the iterative Levenberg-Marquardt 
algorithm, and the parameters in Table 3.2 were recovered to within four decimal 
places. 

Example 2. Two 6-DOF manipulators whose end effectors are rigidly grasping 
together form the 12-DOF closed mechanism in Table 3.4. These parameters were 
used to simulate the movements of this mechanism into 40 different configurations 
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(9j = 0 to 0.5 rad, j = 1 to 6 ). With the initial parameters given in Table 3.3 and 
the simulated joint angles the calibration was performed (also S 3 = 1 fixed). The 
parameters in Table 3.4 were recovered to within four decimal places. 


3.4 Non-Redundant Robot Calibration and Task 
Geometry Estimation 

Next we extend the closed-loop method to situations where the end effector contact 
with the ground has some passive DOFs. Two examples, treated in detail below, 
are a manipulator opening a door (a 1-DOF task) and a manipulator under point 
contact (a 3-DOF task). If t is the number of task DOFs, then the mobility of the 
resultant closed chain is n + t — 6 . For the door-opening task, 6 -DOF non-redundant 
manipulators can therefore be calibrated. For the point contact task, manipulators 
with as few as 4 DOFs may be calibrated. At the same time, the geometry of the 
task is calibrated. 

Since the passive task DOFs are unsensed, they must be eliminated from the 
6 kinematic loop closure equations. Up to 5 unsensed DOFs may be eliminated 
to leave at least one equation for the identification procedure. This elimination is 
simple for the door opening and point contact tasks, but more difficult for arbitrary 
task kinematics. 

3.4.1 Point Contact 

An n-DOF manipulator under point contact is equivalent to grasping a passive spher¬ 
ical ball joint. There are 3 passive DOFs at the endpoint corresponding to orientation 
that are unsensed. Hence the orientation equations in the previous calibration pro¬ 
cedure cannot be used, but the three position equations (3.19) can. Again, define 
the base origin to coincide with the endpoint position. For example, for the 6 -DOF 
manipulator in Figure 3-5, the position of the base coordinates (subscript -1) is coin¬ 
cident with that of the endpoint coordinates (subscript 6 ). The orientation of the -1 
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Figure 3-5: Coordinate description of a manipulator under point contact. 


coordinates is arbitrary with respect to the 0 coordinates, so and a 0 are taken 
as arbitrary constants. Moreover, the orientation of the end effector coordinates is 
arbitrary with respect to the 5 coordinates; hence a 6 (or more generally a n ) is taken 
as an arbitrary constant also. Finally, it is necessary to specify one length parameter; 
for the theorem below, we select a 0 = —1. 

To incorporate just the position equations, we redefine ^ = p‘ = 0 from (3.3) 
and Ax}. = (dx‘, dj/*, dz' c ) from (3.7). The constant parameters , c*o, a n , and a 0 
are removed from <£. Then X, J-, C, and C* are redefined in (3.4)-(3.5) and (3.8) to 
reflect the reduced dimensions. In particular, each column of each Jacobian in (3.12) 
contains only the top vector or first three rows. 





= X i >< b }+! 


(3.27) 


These columns are now interpreted as partial velocity vectors with respect to the 
parameters instead of as screw coordinates [52]. Notice that we have used up three 
kinematic equations in order to eliminate the unmeasured orientation DOFs of the 
endpoint. The identification procedure can then be applied as before. 

As before, the identifiability of the parameters depends on the linear dependence 
of the columns of the Jacobian C*. Because of the form of the columns in (3.27), a 
stronger identifiability condition is derived than (3.15). 

Theorem 4 (Point-Contact Identifiability) The parameters <p are unidentifiable 
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Joint 

s (m) 

a (m) 

a (rad) 

6 °” (rad) 

0 

1.694 

0.837 

3.600* 

0.000* 

1 

1.622 

-0.627 

-0.553 

0.100 

2 

1.000* 

-0.100 

0.930 

0.090 

3 

-0.430 

-0.430 

2.040 

0.100 

4 

0.540 

-0.600 

-0.150 

0.050 

5 

-0.560 

-0.550 

1.350 

0.040 

6 

-1.693 

0.711 

-1.860* 

1.700 


Table 3.5: The initial D-H parameters of a 6-DOF manipulator under point contact. 


Joint 

s(m) 

a (m) 

a (rad) 

i ot} (rad) 

0 

1.594 

0.737 

3.600* 

0.000* 

1 

1.722 

-0.527 

-0.503 

0.120 

2 

1.000* 

0.000 

0.530 

0.000 

3 

-0.330 

-0.330 

2.300 

0.000 

4 

0.440 

-0.440 

-0.900 

0.000 

5 

-0.660 

-0.550 

1.200 

0.000 

6 

-1.793 

0.911 

-1.860* 

1.825 


Table 3.6: The actual/calibrated D-H parameters of a 6-DOF manipulator under 
point contact. 


if and only if there is a constant linear relation among the partial velocity vectors for 
all configurations. That is, there exist constants Cj, k = 1,... ,4, not all zero, such 


that 


n—1 


0 = X c j z j- 1 + X C j z j-1 X b‘ + cjxj + c jXj X b* 


j=0 


3=1 


3+1 


(3.28) 


3=1 


for all i = 1 ,... , to . 


Again, trivial modifications to the theorem can be made if a different length parameter 
than a 0 is fixed. 

We now simulate a 6-DOF manipulator under point contact; the actual parameters 
are given in Table 3.6. The arbitrary constant parameters are marked with a *, and in 
particular we have chosen s 2 = 1 as the fixed length parameter. This entire mechanism 
was simulated in 30 distinct configurations (#i, 0 2 , #3 = 0 to 0.5 rad), starting with 
the initial guesses in Table 3.5. The parameters in Table 3.6 were recovered to within 
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Figure 3-6: The coordinate description of a manipulator connected to a hinge joint. 


four decimal places. 


3.4.2 Opening a Door 

The hinge joints of a door define a rotary axis. Since the endpoint coordinates are 
arbitrary, it is convenient to make z l n coincident with the door’s axis (Figure 3-6). We 
also position the base coordinates at the endpoint coordinates, and let z'_ x coincide 
with z l n . The door hinge angle 9\ measured about z*_ x is unknown, and the orientation 
equation relating to this rotation must be eliminated from the calibration procedure. 

To begin, the position equations are the same as before: p* = 0 from (3.3) and 
Ap* = (dx' c , dy\,dz x c ) from (3.7). The endpoint orientation is given by 

K = R^dx^RyidyDRM) (3.29) 

where dx' c and dy' c are infinitesimal and 0* is finite. Expanding the first column of 
(3.29) and neglecting the second order terms, one finds that 0’ = atan2(R*^ 21 j, 
where the indices denote the elements of the rotation matrix R*. The desired varia¬ 
tions dx\ and dy' c are extracted from R’R z (0’) r = R x (9x*)R v (9y‘). The computed 
endpoint location is then given by the 5-vector x^ = (dx' c , dy' c , dz l c , dx' c , dy' c ). Thus 
one kinematic equation has been used up in order to eliminate the unmeasured door 
hinge angle. 

As before, one length parameter needs to be specified, say ai = —1. Since z , _ 1 
aligns with the door hinge and with z* n , then B°^ and 0 O can be arbitrarily set to 
zero. <£_ is adjusted to eliminate a x and 6 ^*. Finally, X, T, C, and C* are redefined in 
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Joint 

s (m) 

a (m) 

a (rad) 

0°” (rad) 

0 

1.694 

0.837 

3.774 

* 

1 

1.622 

-0.627 

-0.553 

0.080 

2 

1.000 

-0.100 

0.930 

0.090 

3 

-0.430 

-0.430 

2.040 

0.900 

4 

0.540 

-0.600 

-0.150 

0.050 

5 

-0.560 

-0.550 

1.350 

0.040 

6 

-1.693 

0.711 

-1.859 

1.400 


Table 3.7: The initial D-H parameters for a 6-DOF manipulator opening a door. 


Joint 

s (m) 

a (m) 

a (rad) 

9 oit (rad) 


1.594 

0.737 


* 

1 

1.722 

-0.527 


0.000 

2 

1.000 

0.000 


0.000 

3 

-0.330 

-0.330 


0.000 

4 

0.440 

-0.440 


0.000 

5 

-0.660 

-0.550 


0.000 

6 

-1.793 

0.911 

-1.459 

1.825 


Table 3.8: The actual/calibrated D-H parameters of a 6-DOF manipulator opening a 
door. 

(3.4)-(3.5) and (3.8) to reflect the reduced dimensions of a£. The kinematic calibration 
procedure may then be applied. 

Once again, identifiability is related to the rank of C*. Theorem 3 applies to the 
door-opening case, as do the various sources of parameter ambiguity discussed in the 
previous sections. 

Next we simulate a 6-DOF manipulator grasping a door with a hinge joint. The 
D-H parameters are given in Table 3.8; arbitrary parameters are marked by a *, 
including the fixed length s 3 = 1. This entire mechanism was simulated in 40 distinct 
configurations (0\ = 0 to 0.5 rad), starting with the initial guesses in Table 3.7. The 
parameters in Table 3.8 were recovered to within four decimal places. 
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3.4.3 Identifying Arbitrary Task Kinematics 

We now discuss how the algorithm presented for closed-loop kinematic calibration 
readily generalizes to identifying arbitrary task kinematics. As mentioned earlier, the 
chief difficulty is eliminating the unknown DOFs from the environment kinematics. 
This may be achieved by determining the unknown DOFs in terms of the known ones 
(and the kinematic parameters). For instance, in calibrating a system comprised of 
a robot opening a door with a handle, both the door hinge angle and the handle 
angle may be determined in terms of the known manipulator joint angles. Once all 
of the DOFs are determined, the iterative identification algorithm presented above 
is directly applicable to identifying the kinematic parameters. In particular, the 
over-determined system of equations (3.8) may be solved by the iterative Levenberg- 
Marquardt algorithm. 

Determining the unknown DOFs may proceed as follows: (a) using the nominal 
model of the robot, compute the location of the endpoint at a specific configuration, 
then (b) notice that this endpoint also locates the endpoint of the environment kine¬ 
matic chain, and finally (c) using the nominal kinematics of the environment calculate 
the inverse kinematic solution of the endpoint position given in step (a). The result¬ 
ing joint angles are the unknown DOFs. The inverse kinematics of step (c) must 
in general be performed iteratively (for example, with [49]). Notice that a nominal 
model of the kinematics is required. Each iteration of kinematic calibration algorithm 
presented in Section 3.3 improves the nominal model. Thus, the above determination 
of joint angles must be performed anew for each step of the kinematic calibration 
iteration. 


3.5 Discussion 

We have presented a new kinematic calibration method that does not require end¬ 
point measurements or precision points. By forming manipulators into mobile closed 
kinematic chains, we have shown that consistency conditions in the kinematic loop 
closure equations are adequate to calibrate the manipulator from joint angle readings 
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alone. This closed-loop kinematic calibration method is an adaptation of an iterative 
least squares algorithm used in calibrating open chain manipulators. 

Identifiability of the kinematic parameters of the closed loop was reduced to in¬ 
specting the rank of the Jacobian matrix C. Rank degeneracies were in turn studied 
with the screw coordinate interpretation of the columns of the Jacobians C‘. Specif¬ 
ically, a requirement that there be no constant linear relation among the local link 
x*- and z’• axes accounts for all singularities when there are no passive DOFs. Closed 
mechanism with passive DOFs must be studied on a case by case basis for identifia¬ 
bility. 

Three tasks were studied in detail: (1) fixed endpoint, (2) point contact, and (3) 
opening a door. Nevertheless, the method readily generalizes to a large class of robot 
tasks. The main requirement is that there be positive mobility in the closed chain; in 
general, the sum of manipulator DOFs plus the passive DOFs of the endpoint con¬ 
straint must exceed 6. Fixed endpoint constraints generally require redundant arms to 
achieve positive mobility. An equivalent scenario is two manipulators rigidly attached 
at their endpoints with combined DOFs greater than 6; thus two non-redundant arms 
could be calibrated together. When passive endpoint constraints are allowed, single 
non-redundant arms may be calibrated as well; for example, under point contact the 
manipulator only requires 4 DOFs. In principle, up to 5 passive DOFs can be al¬ 
lowed in the endpoint. For every passive DOF, one of the six kinematic loop closure 
equations is used up to eliminate the unknown joint angle; this procedure is akin to 
mechanism synthesis. 

In our method, it is necessary to specify one length parameter to scale the mech¬ 
anism. An independent means for measuring this parameter is required. This is a 
feature of other kinematic calibration methods as well, such as those using theodolites 
|53], 

Another issue with our method is how to handle the forces encountered when mov¬ 
ing the manipulator with an inaccurate kinematic model under endpoint constraints. 
It is beyond the scope of this thesis to present a detailed solution, but an appropriate 
force control procedure must be implemented. The task is made easier if the joint 
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actuation is inherently compliant. For example, one could drive only as many joints 
as there are degrees of mobility in the loop. Alternatively, one could drive all joints 
using the initial guess at the kinematic parameters to calculate the constrained joint 
motion, and then rely on the joint compliance to allow for error in the commanded 
joint trajectories. Though more complicated, this latter method is better, as the for¬ 
mer method could lead to the drive joints jamming at singularities. In Chapter 4 we 
apply the method to calibrating two fingers of the Utah-MIT Dextrous Hand which 
are rigidly attached at the endpoints to form a closed-loop mechanism with 8 DOFs 
(see also [54]); the fingers are moved manually as a simple remedy in this particular 
case. 

Although the joint offsets were the only non-geometric parameters modeled 
in this chapter, in principle more complicated non-geometric joint models could be 
calibrated. For example, in Chapter 4 (see also [54]) an additional scaling factor for 
the joint angle sensors is required and can also be identified. 

Our closed-loop method to kinematic calibration represents a departure from the 
typical dichotomy found in robotics between model building and task performance. 
The removal of this dichotomy may generalize to other problem areas in robotics. In 
this chapter the models of the task and the manipulator are improved while performing 
the task. In Chapter 5 we will show how an uncalibrated stereo vision system can 
be calibrated together with an uncalibrated manipulator (see also [55]). Thus we feel 
that our approach is a step towards true autonomy: no special precalibrated endpoint 
measurement device — or external ‘teacher’ — is needed. 
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Chapter 4 


Closed-Loop Kinematic 
Calibration of the Utah-MIT 
Dextrous Hand 

4.1 Introduction 

The previous chapter proposed a technique for calibrating closed-loop kinematic 
chains formed by dual or redundant manipulators. This technique differs from con¬ 
ventional calibration schemes [34] in that it does not require special endpoint-sensing 
equipment. The present work will experimentally verify this closed-loop calibration 
technique by calibrating the Utah-MIT Dextrous Hand [56] [57]. 

Briefly, the closed-loop kinematic calibration method is described as follows. Con¬ 
sider a finger of the hand system opposing the thumb, such that the fingertips are 
rigidly connected (Figure 4-1). As each finger has four degrees of freedom (DOF) an 
8-DOF closed loop is formed; this loop has in general a mobility of two. Observe 
that fixing only two joint angles uniquely defines the configuration of this mechanism 
(except for the possibility of multiple inverse kinematic solutions). Thus, the other 
six joint angle sensors are redundant. This sensor redundancy may be exploited to 
estimate the kinematic parameters. Specifically, the loop consistency equations for 
a given configuration give three position and three orientation equations containing 
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Figure 4-1: Hand with two fingers opposed, adapted from [4]. 


the unknown kinematic parameters. Moving the fingers into different configurations 
while maintaining the same contact condition provides six additional equations per 
pose. Potentially these equations can be solved for the kinematic parameters, that is 
provided certain identifiability conditions defined in Chapter 3 are met. 

The chapter is organized as follows: (1) outline of calibration procedure, (2) ex¬ 
perimental results, and (3) discussion of future work. Several modifications are made 
to method described in Chapter 3. Initial experiments revealed that the joint angles 
of the hand not only have joint offsets, but also joint scale factors that are difficult to 
determine a priori. For this reason the algorithm is augmented to include the identi¬ 
fication of these joint scale factors. Perhaps the most difficult task in any kinematic 
calibration procedure is determining the initial guess at the kinematic parameters, 
particularly for the base and fingertip transformations. As a partial solution to this 
problem, the more accurately known parameters are fixed while the others are first 
estimated. Finally, singular value decomposition provides a means of dealing with 
parameter ambiguity, and also conveniently produces a measure of parameter observ¬ 
ability. 
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4.2 Method 


4.2.1 Manipulator Kinematics 

Geometric Parameters 

Consider two fingers of the hand connected together rigidly at their tips (i.e., the 
thumb opposed to a finger); see Figure 4-1. This closed kinematic chain has nf — 8 
degrees of freedom. Locate a reference (base) coordinate frame coincident with the 
last joint of the thumb; then number the joints proceeding from that distal joint 
to the palm, and then back out to the tip of the finger, as in Figure 4-1. Let the 
4x4 homogeneous transformation Aj from link j to link (j — 1) be defined by the 
Denavit-Hartenberg (D-H) convention[40] given in Figure 3-2 and symbolically as: 

Aj = Rot(z,9'j)Trans(z,Sj)Trans(x,aj)Rot(x,aj) (4-1) 

where the notation Rot(x, <f>) indicates a rotation about an axis x by (f> and Trans(x, a ) 
indicates a translation along an axis x by a. 

Since the D-H s parameter is not uniquely defined for consecutive parallel axes 
the following Hayati convention[46] (see Figure 4-2) is employed for such axes (i.e., 
for joints 0, 1, 3, 5 and 6): 

Aj = Rot(z,0'j)Trans(x,Sj)Rot(x,aj)Rot(y,aj ) (4-2) 

The position of the last link is computed by a sequence of these homogeneous trans¬ 
formations: 


T c = A 1 A 2 ...A nf (4.3) 

This homogeneous matrix representation of the endpoint is equivalent to the position 
vector (3.1) and rotation matrix (3.2) representation used in Chapter 3. 

The goal is to identify the geometric parameters sj, aj and aj, and also any 
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Figure 4-2: Hayati coordinates, and intermediate x-axis xxj. 
non-geometric parameters that may be included in the kinematic model. 

Non-Geometric Parameters 

The non-geometric effects on the kinematic model potentially include bearing play 
and joint angle sensor error. Although parametric models of both of these processes 
may be included we choose to only model joint angle sensor error. In particular, the 
D-H joint angle is presumed to be related to the sensor reading as 

+ ( 4 - 4 ) 

Thus, we wish to identify the constant joint angle offset and the joint scale factor 
kj. The joint scale gives the calibration factor for the measured electronic signal from 
the Hall effect joint sensors. 

4.2.2 Unidentifiable and Identifiable Parameters 

If we try to identify all the k{ there is a trivial solution ki = 0, which will satisfy any 
set of joint angle data. To avoid this difficulty, fix a joint scale factor that can be 
measured independently. For the hand we choose k? = 1.0. Likewise, one link length 
must be known. This length defines the length scale of the closed mechanism. For 
example, in Section 4.3 we set s 0 = —1.3. 

All of the potentially identifiable kinematic parameters are placed into a single 
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Figure 4-3: Coordinate assignment for loop of two fingers attached at their endpoints, 
vector: 


V = 




(4.5) 


where s = (si,s 2 , ...) r etc. 


4.2.3 Base Coordinate Assignment and Endpoint Location 

The two connected fingertips may be viewed as a single effective endpoint link that 
connects the most distal joint of the finger to the first joint of the kinematic chain 
(i.e., the most distal joint of the thumb); see Figure 4-3. Defined as such, the endpoint 
always has zero orientation and position relative to the base coordinates (i.e., those 
coordinates that were defined to be aligned with the thumb’s most distal joint axis). 
Notice that in addition to the robot kinematic parameters we will also identify the 
D-H parameters of this effective link which completes the loop. 

4.2.4 Endpoint Location Error Calculation 

With a sufficiently good initial parameter estimate the computed endpoint location, 
T c , differs only slightly from zero (the base coordinates). The endpoint error com¬ 
putations can therefore be simplified as follows. Let the computed position (i.e., 
the fourth column of T c ) be represented by (dx Ci dy c , dz c ) T . Similarly, let the calcu¬ 
lated orientation R c (that is, the upper left 3x3 matrix of T c ) be represented by 
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infinitesimal xyz-Euler rotations: 


R x (dx c )Ry(dy c )R z (dz c ) 


1 

-9z c 

dy c 

9z c 

1 

—9x, 

-dy c 

9x c 

1 


(4.6) 


where the right hand side of (4.6) is computed by directly expanding the left hand side 
and ignoring second order differential terms. Thus, the modeled endpoint location, 
evaluated at the ith joint configuration is represented by a six vector: 


£ = ®U£’,£) = {dx c ,dy c ,dz c ,dx c ,dy c ,dz c ) T 


(4.7) 


directly computed from T c . Thus, the explicit form of (3.3) is given by (4.7). 

As stated, the actual position and orientation of the endpoint is taken to be zero, 
thus the endpoint error is given by (3.18). 


4.2.5 Iterative Identification 

Iterative identification may proceed as in Section 3.2.2. 


Differential Kinematics 

At the ith joint configuration £* the first differential of (4.7) is given by a similar 
relation to (3.8): 


A** = ^=r 

96 


A$ + -——As + -Jr—A a + ——Aa + -gr 
9 s 9& 9a 9k 


Ale 


(4.8) 


where again A®* = 0 — ag and As = s — s° etc. By denoting the combined Jacobian: 


& 



9ag 

9g£ 

dag. 

dag 

dag 

M 

9s 

9g, 

9a 

~9k 


(4.9) 
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equation (4.8) may be expressed more concisely as 


Ax* = C'Ay 


(4.10) 


Jacobian Calculation 

From (3.12) the D-H parameter Jacobian elements are: 


and 
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(4.12) 


The Jacobian columns for parameters of the alternate Hayati convention are found 


analogously to be: 
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coli = 


xx j x b’ +1 

dx * 

, col,- = 

y i x b i+1 

XX* 

’ J da 

l- 


(4.14) 


where xxj stands for the local x-axis just prior to the last rotation about the jth 
y-axis by ctj (see Figure 4-2). 

In both parameter conventions the columns of the Jacobian with respect to the 
joint scale factors are: 

P x b* 

z i-i 


1 A i 


(4.15) 


4.2.6 Data Collection and Parameter Estimation 


As in the previous chapter, a series of n configurations of the actual mechanism 
provides n sets of joint angle measurements 0*, and 6n equations of the form (4.10). 
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With the C l defined in (4.9) redefine C as 


C 


C 1 



(4.16) 


A solution may be found iteratively using the parameter error estimates (3.6), and 
<p = <p J0 + A<p. 

As mentioned in Chapter 3, during iteration the matrix C T C may become singular 
at an intermediate parameter set, even though the final parameter set does not have 
a singularity. To avoid the problem the least squares criteria is modified as in Section 
3.2 to be: 


LS' = (A£-CA£) t (A*.-CA£) + A(A£) t (A£) (4.17) 

In contrast to Chapter 3, this criteria is minimized here by using the singular value 
decomposition (SVD) of C, zeroing singular values that are less that p percent of the 
maximum singular value, and then implementing the generalized inverse from the 
SVD matrices [58]. The value of p implicitly gives A, and it is set high (e.g., p = 5 
percent) initially and reduced once convergence occurs. 

If the kinematics are over-parameterized for the collected joint angle data then 
C will always have zero singular values. The number of near zero singular values 
indicates the number of non-independent parameters. 

Two other variations of the original algorithm are as follows. The poorly known 
parameters (base and tip transformations) are first estimated by assuming that the 
finger parameters are correct. Then all parameters are allowed to freely vary, giving 
the final results. Also, it is often found useful to check that the present update to the 
parameters improves the endpoint positioning error. If it does not then the parameter 
update is repeatedly halved, until the endpoint error improves (see [58]). 
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Joint 

9°” (rad) 

* (in) 

a (in) 

a (rad) 

k 

0* 

0.0000 - 

-1.3000 

0.0000 

0.0000 

1.0000 

1* 

0.0000 

-1.7000 

0.0000 

0.0000 

1.0000 

2 

0.0000 

0.0000 

-0.4500 

1.5708 

1.0000 

3* 

-1.5708- 

2.1200- 

0.4094- 

0.2080- 

1.0000 

4 

0.0000- 

2.2500- 

0.6000- 

-1.5708 

1.0000 

5* 

0.0000- 

1.7000 

0.0000 

0.0000 

1.0000 

6* 

0.0000 

1.3000 

0.0000 

0.0000 

1.0000 

7 

3.1415- 

-0.5000- 

-1.0000- 

-0.2618- 

1.0000 


Table 4.1: 8 DOF mechanism. Initial parameters. A * indicates that the parameters 
are those of the Hayati convention (with the respective units). The parameters marked 
with a — are difficult to measure and were either guessed (for link 7) or roughly 
calculated from known specifications. 

4.3 Experimental Results 

Finger three of the hand was opposed to the thumb (finger 0) by rigidly connecting 
the fingertips (screwed into a common aluminum plate). The nominal kinematic 
parameters for this 8 DOF closed kinematic chain were taken from [57] and are shown 
in Table 4.1. The geometric parameters (i.e., the s, a and a parameters) marked with 
a — are difficult to measure independently. The joint angle offset and scale parameters 
are also not well known. Those geometric parameters not marked with a — are likely 
to be more accurate than parameters identified by any identification scheme relying 
on the limited accuracy of the joint angle sensors. 

A series of 200 configurations of the finger/thumb mechanism provided input joint 
angles for the above identification algorithm. These joint angles are plotted in Figure 
4-4. The joint angles for the thumb are negated so that the identified joint axes are 
in the same direction as defined in [57]. Considerable care was taken to make sure 
that all joints moved as much as possible. 

The calibration was performed with the initial parameters given in Table 4.1, 
the recorded joint angles, and the two fixed parameters 5 0 = —1.3 and kr — 1.0. 
Convergence of the algorithm required p = 0.5 in the SVD based pseudo-inverse, and 
produced the parameters in Table 4.2. 
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Joint 0 off (rad) s (in) a (in) a (rad) k 


0* 

0.2380 

-1.3000 

0.0069 

0.0668 

0.9508 

1* 

0.0850 

-1.5085 

-0.1515 

-0.0480 

0.8596 

2 

-0.0943 

0.1604 

-0.3097 

1.6207 

0.9336 

3* 

-1.7824 

1.8050 

0.5755 

0.2890 

0.9579 

4 

0.1389 

1.7151 

1.3356 

-1.5836 

1.0125 

5* 

0.5464 

1.6307 

0.0292 

0.0737 

0.8577 

6* 

-0.0425 

1.1474 

0.1745 

0.1260 

0.9350 

7 

3.7552 

-0.8322 

-1.4377 

-0.5631 

1.0000 


Table 4.2: 8 DOF mechanism. Calibrated parameters. A * indicates that the param¬ 
eters are those of the Hayati convention (with the respective units). 
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4.3.1 Endpoint Errors 

The calibration may be assessed by measuring the error in closure between the fin¬ 
gertips. This error is simply given by the position and orientation vector Ax* for 
each pose i. Averaging over all poses, the pre-calibration root mean square (RMS) 
position error was 0.5780 in and the RMS orientation error was 0.4645 rad. After 
calibration these errors are reduced to 0.0201 in and 0.0290 rad respectively. 

4.3.2 Geometric Parameter Errors 

As stated, the geometric parameters not marked with a ~ are likely to be more 
accurate than the parameters identified by any identification scheme relying on the 
limited accuracy of the joint angle sensors. These parameter values may be used 
to check the validity of the calibration algorithm. Comparing these parameters in 
Tables 4.1 and 4.2 reasonable agreement is seen. The existing miss-match could be 
partly due to unmodeled kinematics such as joint wobble, and other factors caused 
by machining imprecision. Also, as discussed below, certain parameters may not be 
uniquely determined from the limited accuracy joint sensor readings. 

4.3.3 Non-Geometric Parameter Errors 

It is seen from Table 4.2 that the joint angle offsets and scale factors are an impor¬ 
tant source of error in the kinematic model. The joint offsets marked with a ~ are 
particularly hard to measure, and indeed show the greatest error. The non-geometric 
parameters that are more accessible (i.e., and kj for j — 0,1,6) were carefully 
re-measured after calibration. This post hoc measurement indicated that these non¬ 
geometric parameters were identified more accurately than the initial guesses in Table 
4.1. 
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4.4 Parameter Identifiability 

The ratio of the minimum singular value to the maximum singular value of C provides 
an index of parameter identifiability. This ratio is 0.003, indicating that C T C is not 
singular, but is not well conditioned. Thus, all the parameters are theoretically 
identifiable, but some may be sensitive to measurement noise. As stated, to obtain 
convergence all the singular values less than 0.5 percent of the maximum singular 
value were zeroed. This translated into zeroing seven singular values. Thus, for this 
mechanism, and for this particular joint angle data set, there are seven parameters 
that are close to linearly dependent upon the other 31 parameters. 

It is difficult to determine the geometrical significance of these interdependent pa¬ 
rameters. In theory, linear relations between the instantaneous parameter variations 
may be found from the null space of the Jacobian (the null space is determined by the 
span of the columns of the matrix V that correspond to (near) zero singular values, 
where V is part of the SVD of the Jacobian: C = UDV T ). It is not clear though 
how to interpret these first order relations to determine the source of the parameter 
ambiguities. 

Certain parameter dependencies do have a simple explanation though. For ex¬ 
ample, consider the parameters associated with the most proximal joint of the finger 
and the most proximal joint of the thumb (joints angles 3 and 4). These two almost 
parallel joints are partially decoupled from the rest, in the sense that is possible to 
move them without large movements of the other joints and visa versa. From the data 
(Figure 4-4) it is seen that joints 2 and 3 move the least (only 0.4 rad). Roughly, 
these two joints, being perpendicular to the others, move the plane of the rest of the 
finger/thumb complex. The location of this plane is given by the parameters a 2 and 
a 4 , which are measured along the common normals X2 and X4. Unfortunately, these 
two common normals stay almost parallel for all the data collected (mostly because 
of the small movement ranges of the joints 2 and 3). In the extreme case where they 
are fixed to being exactly parallel then the D-H distances a 2 and a 4 along these com¬ 
mon normals may vary arbitrarily, as long as the difference a 2 — a 4 is kept the same. 


65 



Joint 

S ott (rad) 

* (in) 

a (in) 

a (rad) 

k 

- 1 * 

0.0000 

1.8050 

0.5755 

0.2890 

N/A 

0 

0.1389 

1.7151 

1.3356 

-1.5836 

1.0125 

1 * 

0.5464 

1.6307 

0.0292 

0.0737 

0.8577 

2 * 

-0.0425 

1.1474 

0.1745 

0.1260 

0.9350 

3 

? 

? 

? 

? 

1.0000 


Table 4.3: Identified finger parameters relative to a palm reference frame, where a * 
indicates that the parameters are those of the Hayati convention (with the respective 
units), and a ? means an unknown tip transformation. 

Though these common normals are not exactly parallel to one another, we should 
predict a sensitivity problem in separately identifying a 2 and a 4 from finite precision 
measurements. See Chapter 3 for a more general discussion of identifiability. 

The identifiability of these parameters may be improved by relaxing the endpoint 
constraint to be a point contact (with the endpoint free to orient arbitrarily) instead 
of a rigid contact. In the case just discussed the axes x 2 and x 4 would then no 
longer be constrained to be nearly parallel. The use of a point contact constraint 
(or a perhaps more elaborate rolling contact constraint) has not yet been attempted 
because, as yet, we do not have tactile sensors that could be used to assure the point 
contact. We could build a ball joint to connect the fingertips, but this would not be 
as general, or natural way of proceeding. 


4.5 Common Palm Reference Frame Conversion 

Once a finger is calibrated against the thumb, as just described, it is necessary to 
convert the identified parameters into a reference frame located on the stationary 
palm of the hand, similar to [57]. This reference frame is the frame labeled 2' in 
Figure 4-3, that is, the frame found after the rotation of local frame 2 by the joint 
angle 0 3 . The finger’s parameters relative to this palm reference are therefore given 
by the parameters labeled 3 and greater (except 0 3 ), and are re-written in Table 4.3. 
The tip transformation is not identified by the above method. It could be identified 
if a point contact constraint were used instead of a fixed constraint. For practical 
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Joint 

(rad) 

* (in) 

a (in) 

a (rad) 

k 

0 

1.7824 

0.0000 

0.3097 

-1.6207 

0.9336 

1 ’ 

0.0000 

-0.1604 

0.0000 

0.0000 

N/A 

1 * 

0.0943 

1.5085 

0.1515 

0.0480 

0.8596 

2 * 

-0.0850 

1.3000 

-0.0069 

-0.0668 

0.9508 

3 

? 

? 

? 

? 

1.0000 


Table 4.4: Identified thumb parameters relative to a palm reference frame, where a * 
indicates that the parameters are those of the Hayati convention (with the respective 
units), and a ? means an unknown tip transformation. 

purposes the tip D-H parameters are taken to be 0°** — 0, s = 0, a = 0.735 and 
a = 0 [57]. Finally, notice that the parameter set labeled —1 locates the first axis of 
the finger relative to the thumb. 

The thumb parameters are found by following the kinematic chain from the palm 
frame 2' backwards to the ‘base’ (actually, the fingertips). This simply entails re¬ 
versing the sign on the respective 0 , s, a and a parameters. Also, in the case of the 
Hayati convention parameters the order of transformations must be changed to: 

Aj = Rot(z,0'j)Rot(y,aj)Rot(x,a,j)Trans(x,Sj ) (4.18) 

The thumb parameters are thus calculated to be those shown in Table 4.4. Notice 
that an additional s translation, s^, is required to provide the translation of —0.1604 
along the thumb joint 1 axis. 

When the whole calibration procedure is repeated with other fingers opposing the 
thumb then the identified palm reference frame 2’ may differ from finger to finger. 
That is, this frame may be arbitrarily displaced or rotated on the thumb’s first joint 
axis (labeled z 2 in Figures 4). One such reference frame may be defined as the 
‘common reference frame’ and the others related to it. The unknown translation and 
rotation on z 2 (between any two palm reference frames) may be estimated by locating 
the thumb in a fixed position and calculating its endpoint in both reference frames. 
The difference between these endpoint locations gives the two unknown quantities. 
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4.6 Discussion 


Studies investigating the kinematic calibration of the Utah-MIT Dextrous Hand have 
been presented. The approach was to use the closed-loop kinematic calibration tech¬ 
nique developed in the previous chapter. Several modifications of the original algo¬ 
rithm were necessary. First, it was found that joint angle scale factors had to be 
included as parameters to be identified. Second, it proved to be important to iden¬ 
tify initially the endpoint and palm (base) transformations by assuming the initial 
guesses at the finger kinematic parameters were correct. This provided improved ini¬ 
tial estimates of these parameters, which were otherwise tricky to measure. Once this 
first step was completed the full scale non-linear parameter calibration technique was 
used. Finally, SVD provided a convenient way of resolving the parameter ambiguity. 

The experimental results indicate that the kinematics of the finger/thumb com¬ 
plex can be identified by the proposed closed-loop kinematic calibration method. 
Endpoint positioning error was improved by over an order of magnitude. Further, 
the identified values of the parameters that were also accurately known a priori were 
in close agreement with these a priori values. 

Certain parameters were not uniquely determined from the collected joint angle 
data. This parameter ambiguity, while clearly not a problem for the particular finger 
contact situation studied (see RMS error), may cause troubles when the identified 
kinematics is used for very different finger configurations. The parameter ambiguity 
is principally due to the lack of joint movement allowed by the fixed constraint (it 
may also be due to limited joint sensor resolution, and other unmodeled kinematics). 
In the future we intend to use a point contact constraint (as described in Chapter 3 
and [39]). This will produce an 11 DOF mechanism with higher mobility. The natural 
way to implement this is to use tactile sensors to assure that the fingers stay at a 
point contact, while allowing arbitrary orientation of the fingertips. This experiment 
awaits the completion of tactile sensors for the fingers. Other more complex endpoint 
constraints can also be explored. In fact, the above technique may be extended to 
identifying the geometry of an object grasped between the fingers (see Chapter 3 and 
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Chapter 5 


Autonomous Calibration for 
Hand-Eye Coordination 


5.1 Introduction 

Hand-eye coordination is a particularly demanding task because it requires the con¬ 
sistency of two separate sensing systems — the robot manipulator and the stereo 
vision system. It is the intention in this chapter to address the issue of how these two 
systems may be calibrated autonomously. By calibration it is meant the determina¬ 
tion of the parameters of the internal models of both the camera and arm systems. 
As in the the previous chapters, the emphasis is on autonomy. Thus, only the robot’s 
internal joint angle and camera image sensors are permitted for the calibration. 

5.1.1 Motivation 

Autonomy is a particularly important property for a robot that must function outside 
of controlled laboratory conditions. It is inevitable that a robot will have its base 
moved, links bent, cameras misaligned or be otherwise damaged. In such situations 
it would be desirable not to have to resort to the use of special purpose calibration 
equipment to update the model used for robot control. In fact, an ultimate goal 
would be for the robot to be able to calibrate its internal model in real time. 
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While in certain engineering applications the goal of autonomy may be sacrificed 
in favor of simplicity, it was pointed out in Chapter 2 that humans have no such 
choice in calibrating their sensorimotor system. Thus, a second motivation for study¬ 
ing autonomous calibration derives from an interest in understanding the human 
sensorimotor system. 

5.1.2 Previous Research 

In the domain of robot dynamics autonomous calibration has essentially been achieved, 
although the kinematics must be assumed to be known (for example [16]). In par¬ 
ticular, it is possible to estimate the inertial parameters that define the various links 
by only using internal joint torque (current), position, and velocity measurements. 
This idea has actually been made to operate in real time model-based adaptive con¬ 
trol schemes[59]. The success of inertial estimation is based on the observation that 
the suitable combinations of the inertial parameters enter linearly into the dynamic 
equations. 

In contrast, autonomously determining the static relationship between the internal 
joint angle sensors and the manipulator endpoint position — the kinematic model 
- has not been as successful as autonomous dynamic estimation (see review [34]). 
Typically, researchers have viewed the manipulator as a positioning device — that is, 
an open-loop kinematic chain. This view demands that the endpoint be measured in 
addition to the joint angles to infer the kinematic parameters. Therefore, autonomous 
calibration is not possible. 

If instead the manipulator is viewed as a device to interact with the environ¬ 
ment then autonomous calibration is possible, as described in Chapter 3 (see also 
[37] [39] [38]). The basic observation is that the manipulator may form a mobile closed- 
loop kinematic chain when performing a task. The internal model structure, the 
knowledge of the type of task constraint, and the internal joint angle measurements 
collected while in a number of configurations provide enough consistency equations to 
solve for the kinematic parameters. This technique of closed-loop kinematic calibra¬ 
tion is quite general. As an instance that is most relevant to the hand-eye calibration 
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problem, we recall from Chapter 3 that two uncalibrated manipulators may calibrate 
one another by moving while rigidly grasping (a O-DOF task) a common object. 

Three observations are worth stressing: 

• The knowledge of the task constraint (e.g., two robots gripping together) re¬ 
places the need for an external sensor. 

• The redundancy of the sensing systems (e.g., two arms) with respect to the task 
enables the various redundant components to move while performing the same 
task. 

• The a priori model structure knowledge allows one to form a number of con¬ 
sistency equations (the kinematics) which may be solved for the kinematic pa¬ 
rameters. 

These three observations serve as a basis for extending autonomous manipulator 
calibration to complete hand-eye calibration. First we review relevant vision system 
calibration techniques. 

The conventional methodology for camera calibration is to move a number of 
known precision points into the field of view of the cameras and infer the camera 
calibration from the given points in space (see review [60]). 

One effective approach is to form a look-up table from known rays (obtained from 
two planes of points in space) to recorded image locations [61] [62], and then use 
splines to do local interpolation. Look-up table approaches need external calibration 
points, and thus they must be disregarded for autonomous camera calibration. 

Various model based approaches have been used for camera calibration. In com¬ 
puter vision and graphics the pinhole camera model has been used extensively [63] 
[64] [65] [66]. This model may be augmented to account for lens distortions [67] [68] 
[69]. Since the pinhole camera model is non-linear in its parameters there have been 
various proposals to make the calibration equations linear [31] [70] [69]. These meth¬ 
ods are important because they provide initial guesses at the parameters that general 
non-linear optimization methods require. Other empirical polynomial interpolation 
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models have also be used [61] [71] [62], their only advantage being that the parameters 
enter linearly into the equations. 

While the above camera modeling techniques and parameter identification meth¬ 
ods are relevant to autonomous calibration, most of this work is predicated on the 
assumption that there are external calibration points available. There are a few no¬ 
table exceptions. For example, early photogrammetric engineering work [72] and more 
recent robotics research [73] have demonstrated that the camera parameters may be 
recovered by moving the cameras while viewing arbitrary unknown points in space. 

Finally, there is a considerable body of literature that addresses the problem of 
registering the calibrated vision system’s coordinates with respect to the robot base 
(see review [34]). Especially interesting is the work in [74]. Their technique is to 
determine the camera to hand transformation simultaneously to the robot parameters 
by viewing an arbitrary fixed point in space. The internal camera parameters are 
calibrated beforehand by viewing a precision calibration jig. 

5.1.3 Towards Autonomous Hand-Eye Calibration 

As can be seen the calibration of a robot manipulator/vision system is typically 
based on a ‘divide and conquer’ principle. None of these approaches may be made 
autonomous. Further, there is no guarantee that once the various separately cali¬ 
brated components are assembled they will be consistent. This is a very important 
point since we know that models and sensor readings are inaccurate. We thus pro¬ 
pose that a solution is to calibrate the models of the manipulator and two cameras 
simultaneously while performing the task of interest — hand-eye coordination. This 
may be done as follows. 

Recalling the three observations made concerning closed-loop manipulator cali¬ 
bration, it is remarked that hand-eye calibration fits into this framework. First, a 
manipulator and a vision system may sense the location of the same point in space; 
and thus, the total robot sensing is redundant. Second, if the task is defined as the 
cameras tracking the hand then there is a closed kinematic loop formed. This task 
constraint replaces the need for external calibration points. Finally, because we may 


73 



assume a priori knowledge of both the camera and arm kinematic model structure 
it is possible to write out consistency equations of the closed loop in a number of 
configurations, and thus solve for the parameters. We will develop this idea in what 
follows. 


5.1.4 Outline 

As stated, the purpose of this chapter is to extend the closed-loop calibration approach 
developed for calibrating robot arms to calibrating a complete robot — with a stereo 
vision system in addition to the manipulator. The stereo system is assumed to have 
one axis of rotation per camera, but is otherwise taken to have an arbitrary geometry. 

An uncalibrated stereo camera system will be made to track a point on the hand 
of an uncalibrated arm. There are at least two distinct approaches to forming the 
closed-loop calibration equations. The first is to formulate a model for the manipu¬ 
lator relative to each camera separately, and measure the position error in 2-D image 
plane coordinates. The calibration would proceed by collecting data from manip¬ 
ulator/camera movements and minimizing the image plane error in both cameras. 
The second approach is to directly model the position of the end effector given by 
the stereo calculation (from the image data). The calibration can then proceed by 
minimizing the end effector error between the manipulator and stereo models. Of 
these two approaches the second is chosen because it seems more natural to minimize 
the task space error. In addition, using the second approach enables one to formu¬ 
late the iterative identification equations more simply; in particular, the manipulator 
Jacobians developed in Chapter 3 are directly applicable. 

It is assumed that the point that is to be tracked may be unambiguously located 
on both camera images. This generally non-trivial correspondence problem [66] may 
be solved here because of two additional constraints. First, it is known that the hand 
is moving relative to the background; therefore, it is possible to disambiguate the 
hand image from the background. Second, a convenient point that can always be 
unambiguously located on the hand may be used (e.g., the tip of a pointer). 
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5.2 Model Definitions 


5.2.1 Manipulator Model 

Consider an arbitrary manipulator with nj degrees of freedom. Let the 4x4 ho¬ 
mogeneous transformation Aj from link j to link (j — 1) be defined by the Denavit- 
Hartenberg (D-H) conventional)] given in Figure 3-2 and (4.1). 

For convenience we define the base of the manipulator to coincide with a head 
referenced coordinate system that is coincident with the left camera’s axis of rotation 
(see next section). The position of the last link is related to these base coordinates 
by a sequence of D-H transformations defining the kinematic model: 

T c = A 0 A 1 A 2 ...A nf ( 5 * 1 ) 

Since the model must only locate a point on the end effector the last axis skew is 
unnecessary, i.e., let a n/ = 0. 

The position of the point P to be tracked is thus given by the model: 


x 

1 


Aq A\ Ai ... An^ Q. 


(5.2) 


where the • emphasizes that it is the position of P modeled by the manipulator system, 
and o = (0,0,0,1) is the location of the point P in the hand based reference frame. 
All of the unknown manipulator kinematic parameters are placed into an array: 


T\T 


l = &„,* T ,a T ,a T ) 


where s — (s x , s 2 » — ) T etc. 
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5.2.2 Visual System Model 

The model of a stereo camera system may be decoupled into a purely geometric part 
giving the relative orientation and position of the cameras, and a part modeling each 
camera’s projection of points in space. The parameters for these two portions are 
respectively called the external and internal camera parameters. 

The internal camera model used is the standard pinhole camera model [66] [75]. 
More refined parametric models including lens distortions [67] [68] [69] may also be 
incorporated without changing the general approach. Let the effective focal point 
distance be denoted /, and the projected point P in the image plane be given by the 
pair (u,v). Further define local camera coordinates to have x and y axes parallel to 
the camera ( u , v) grid and have an origin at the focal point (see Figure 5-1). Thus, the 
coordinates (z R ,yj j, zr) of a point P expressed in the right camera’s local coordinate 
system is given by the standard projection equations: 

W(~/r) = x R /(u R + u°J f ) = y R /(v R + v% f ) (5.3) 

Notice that provision is included for unknown offsets (ur* ,v R *) between the image 
plane readings and the optical center of the camera. Analogously, the left camera 
equations are: 


W(~/i) - X L /(u L + U° L ff ) = VL/(vL + vl ff ) (5.4) 

The external geometric model of the two camera system can be represented by 
D-H transformations. To distinguish the camera D-H parameters a tilde (e.g., So) is 
used. We assume that each camera has one degree of rotation about a fixed axis with 
a joint angle sensor (9r = 9 2 and Or = 9\ for the left and right cameras respectively) 
as defined below. It is convenient to start the kinematic chain at the right camera’s 
local coordinate system — that is, the frame located at the right camera’s focal point 
and having its x and y axes parallel to the right image plane («, v) coordinate grid. As 
mentioned in the previous section, the ‘base’ coordinates are assumed to be located at 
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the left camera axis of rotation. Thus, the transformation of a point x = ( x,y,z , l) r 
in base coordinates to the local coordinates of the right camera x R = (xr, y R , zr, l) r 
is given by (see Figure 5-2): 


Xr = A 0 A 1 X = Tr'x (5.5) 

where the only variable parameter is §[ = 0 R + , the right camera rotation. We 

also define the opposite transform as x = TrXr. The analogous left camera coordinate 
system may be located by a further two D-H transformations (see Figure 5-3): 


x = A 2 A 3 Xl = T l Xl 


(5.6) 


where the only variable parameter is 0' 2 = 0^ + 6° 2 *, the left camera rotation. The 
parameters S 3 and translate and orient about the left optical axis so that the left 
coordinate system lines up with the ( u, v) image grid, and is located at the focal point. 
The parameters 03 and &3 are redundant and are taken to be fixed zero quantities. We 
also define Rr and Rr as the upper 3x3 rotation matrices of Tr and Tr respectively. 

To solve the stereo camera equations it is convenient to define the vector p from 
the left focal point to the right focal point, the vector 1 from the base frame origin to 
the left focal point, and the two internal parameter vectors as follows: 
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1 

1 




-RlPl , 
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T L 0 
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UL + U 0 L f 
vl + v° L fi 
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Figure 5-1: The internal camera model. 
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where again o = (0,0,0,l) r . Notice that U£, and u# are respectively vectors from 
the left and right focal points along the line of sight. Thus, the point P in base 
coordinates is simply: 


X = CU£ + 1 (5.11) 

where the • emphasizes that it is the position of P modeled by the camera system. 
The scalar c is given by: 


(5.12) 

The linear 
All of the 

i = (£ff,l T ,a T ,« T ,i T ) T 

where i = ( f R , f L , u°J f , v °^ f , u° L ff , vl ff ) T and S = (Si, S 2 , ...) r , etc. 


c = i!** u *- u * x P = r(p) 

Ufl X U L . Ur X U £ 


where ‘x’ denotes vector cross product, and *•’ denotes inner product, 
operator T(-) has been defined here as it will be useful in Section 5.1.1. 
unknown camera parameters are placed into an array: 
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Figure 5-2: The right camera axes to base D-H transformations. The first transforma¬ 
tion locates the line of action of the right camera rotation. The second transformation 
locates the base coordinate z-axis, which is also the left camera axis of rotation. 



Figure 5-3: A stereo camera system attached to a manipulator. L and R indicate the 
left and right cameras. The left and right cameras rotate about Zi and z~o respectively. 
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5.2.3 The Closed-Loop Model 

Comparing equations (5.11) and (5.2) it is seen that the difference defines the closed- 
loop kinematic equations of the hand-eye system: 

0 = i(<£) — = A®(^) (5.13) 

->T -T T 

where we have defined ^ = (<£ ,<£ ) as a concatenation of all of the parameters 
to be identified. Note that only unknown parameters need be included in <f>. For 
example, if all the camera parameters were known and not included in <f> then the 
calibration method described would be a standard manipulator calibration scheme, 
with endpoint visual sensing. 

Also, recall that the base coordinates of the manipulator were defined to corre¬ 
spond to the local coordinates of the left camera rotation. Thus, we see that the 
manipulator base coordinates are the camera coordinates with axes x 1} y 1? and z 1 . 


5.3 Model Identification Procedure 

As the cameras track the point P at discrete locations the joint angle and image 
plane sensory information is recorded. The data recorded at the ith configuration are 
placed into a single array: 




(5.14) 


At the ith configuration of the hand, one vector equation of the form (5.13) may be 
written 


0 = A x* = A *(«*,$ (5.15) 

where in additional to the functional dependence on the model parameters the de¬ 
pendence upon the ith data array u* is explicitly shown. As a short form Ax' will be 
used, where the functional dependence on the ith data array u l will be understood 
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by the super-script i. After moving the hand and cameras into m distinct locations 
3 m scalar equations are generated. These equations must be solved in order to find 
the optimal parameter set <£, completing the calibration. 

5.3.1 Iterative Identification Technique 

We use the same iterative method introduced in Chapter 3 to search for the solution 
By expanding (5.13) with a Taylor’s series about an initial guess and neglecting 
higher order terms, the following linearized form is obtained: 


Ax' 


d± 



(5.16) 


where A<£ = 4^ — $• The Jacobian C ’ may also be written as: 


d& 

■ dss.' 

,dx l 

■ <9x* | dx' 
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~di 
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'~da 

~da ~d[ 

ds 

'da 
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1 di 


Recall that Ax* is the difference between the location of the point P given by the 
camera model and the manipulator model (computed using 4>q). Until the system is 
calibrated this difference is non-zero. 

Place the equations for each of the loop configurations into one array: 


AX = 


Ax 1 


C 1 

Ax 2 

_ 

c 2 

Ax m 


c m 


A^> = CAfi 


(5.17) 


The Levenberg-Marquardt style non-linear search described in Section 3.2.3 is used 
to iteratively estimate the parameters. To reiterate, the solution which minimizes the 
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modified least squares criteria (AX — CA<£) T (AX — CA<£) + AA<£ T A<£ is (3.17): 


At = (^C + Xiy^AX 


(5.18) 


The current guess at the parameters is iteratively updated using (5.18) and <f> = 

Jacobian Calculation 

The columns of the manipulator Jacobian were derived in (3.27) to be: 


, dx' 
col ^ = 
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(5.19) 


(5.20) 


The partial derivatives of the camera model may be obtained by similar methods. 
First concentrate on the length parameters. Notice that from (3.1) the vectors p and 
1 in equations (5.7) and (5.8) may also be written as: 


p = 3jZj_i + a,jXj 

3=0 
3 

1 = X 3jZj_ i + ajXj 

j=2 


(5.21) 

(5.22) 


where z and x are local camera z and x axes. Thus, the camera model of point P (5.11) 
can be re-written to explicitly show the linear dependence of the length parameters: 


= - £ 5 i[ r ( 2 i-i)u£] + «i[r(xi)u L ] 

3=0 

3 

+ X 1 - r (*i-iK] + a i [x,- - r(x i )u £ ] 

3 =2 


(5.23) 


where we have taken advantage of the linearity of the operator T(-) defined in (5.12). 
It is thus apparent that the column of the Jacobian with respect to a particular length 
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parameter is given by the term in square brackets that it multiplies in the previous 
equation. For example, the columns with respect to Si and S 2 are respectively: 


col. 


dx* 

dsi 


r(z 0 )u L 


dx' 

’ cc ’ li w 2 = 


z"i - r(ai)u L 


(5.24) 


The other parameters may be treated by direct differentiation of (5.11). For 
example, consider the movement of the point P due to a variation in 0 : 2 : 


= Aa 2 [-^-u L + cx 2 x bg] 
0OC2 


(5.25) 


where = U&. The term in square brackets in (5.25) is the column of the Jacobian 
with respect to d 2 . The evaluation of dc/da 2 is straightforward but messy. It involves 
only the partial derivatives of Ul and p which are respectively x 2 x b^ and x 2 x p. 


5.4 Model Identifiability 

Inspecting the form of C defined in equation (5.17) it is seen that the columns of 
C will be independent, and thus the model identifiable (as defined in Chapter 3) if 
and only if there does not exists some fixed linear relation among the columns of the 
Jacobians <7*. Using (5.19), (5.20), (5.24) and (5.25) the following condition condition 
is obtained: 

Identifiability condition: identifiability is guaranteed by checking that 
there be no constant (for alii = l..m configurations) linear relation among 
the manipulator local link x*- axes and z’- axes, their moment vectors 
x*- x b’- and z*-^ x b*-, and the camera Jacobian vectors given in the 
previous section. 

The identifiability of the parameters manipulators was studied in detail in previous 
chapters; thus, we will restrict ourselves to camera related problem situations. 

As a first example of unidentifiable parameters, consider the situation when the 
two camera rotation axes are parallel: zj, = z\ = z*. One might suspect a problem 
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because it is well known that the D-H s parameter is poorly defined for serial link 
kinematic chains with consecutive parallel axes [46]. The columns of the Jacobian 
(5.24) with respect to and J 2 are r(z’)u^ and z* — T(z*)u^ respectively. These two 
vectors are not linearly dependent, so the problem is not as simple as for serial link 
chains. Notice though that z\ = z’ is the base coordinate’s z-axis of the manipulator, 
and the Jacobian vector with respect to the manipulator parameter so is z. Thus, 
there exists a linear relation among these three Jacobian vectors: 

0 = -[z*] + [z* - r(z<)u*J + [r(z*X] Vi (5.26) 

and by the identifiability condition the parameters s 0 , «s x and s 0 are not identifiable 
alone. This becomes a practical problem because it also means that C T C will be sin¬ 
gular in the iteration algorithm (5.18). The solution is to use an alternate coordinate 
convention for the transformation A x . The Hayati convention [46] developed for ma¬ 
nipulators provides such a four parameter system (see Chapter 4). This convention 
cannot be used for all joints because it too has a similar ambiguous parameter when 
there are two consecutive perpendicular axes. 

As a second instance of unidentifiable parameters, notice that the closed-loop 
equations (5.13) of the calibrated model may be written to show explicitly the linear 
dependence of the length parameters: 

0 = )u L ]+ a j [r(x i )u i/ ] 

i=o 

3 

+ J2 s A i j-i - r (zi-i) u i] + - r(x,-)u L ] (5.27) 

j=2 

tlf 

~Y, s A z j-i] + «i[ x i] 

3=0 

where in addition to equations (5.13) and (5.24) the fact that the manipulator length 
parameters enter linearly into the kinematics has been used. All of the terms in square 
brackets are columns of the Jacobian (7*, and thus, the closed loop is unidentifiable by 
the above ‘identifiability condition’. Also from (5.27) it is seen that the source of the 
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trouble is that the loop equations may be scaled arbitrarily and still satisfy the joint 
angle and image plane data. The solution to this problem is to fix one link length 
arbitrarily. This length determines the units of length. If conventional units, such 
as a meter, are desired the calibrated vision system merely has to view the desired 
meter stick and calculate the correction scale factor to its internal units of length. 

As a final example of unidentifiable parameters consider the case when the data is 
not ‘persistently exciting’. For instance, if the hand point P is always to move such 
that it stayed in a plane defined by the right camera axis z' 0 , and this plane happened 
to be coincident with both focal points, then (u# x u^) x z l Q = 0 for all i configurations. 
Thus, r(z* 0 ) = 0 and the column of the Jacobian with respect to Si is identically zero. 
The parameter Si is unidentifiable in this situation. Though this scenario is unlikely 
it does point out the importance of the choice of the configuration data used for 
calibration. Although perhaps only simulation may determine the optimal data set, 
it is possible to study the sensitivity of particular parameters [76]. 


5.5 An Example and Simulation 

In order to clarify the general calibration procedure an example is now presented and 
solved by simulation. Consider the planar 2-DOF manipulator connected to a head 
mounted stereo system in Figure 5-4. In total there are six degrees of freedom in 
the system: two manipulator joints, two camera rotations, and two one-dimensional 
images. The kinematic parameters of interest are in the following order placed into 
the array <£: the three link lengths, three joint angle offsets, three length parameters 
providing the displacement between cameras coordinates, two camera rotation offsets, 
and two focal lengths: 0 = (ao.oijOa,^,^ ,s o ,a u s z ,0l if ,0 2 ff , f L , f R ) T . A 

simulation of this hand-eye system was performed by using the ‘actual’ parameter 
values 

0 = (1-0,1.0,1.0, —.8,0,0, .1, .36, .1,0,0, .05, .05) r (lengths are in meters and angles in 
radians). Joint angle data was generated by moving the four rotational joints 0 X , 0 2 , 
9 lf and 0 2 over a trajectory starting at 6 X = 0, 0 2 = 0, 9\ = 0, and 0 2 = 0, and covering 
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Figure 5-4: Example hand-eye system used for simulation purposes. 

a joint space volume defined by increasing each joint three times in increments of 0.1 
rad- i.e., 3x3x3x3 = 81 distinct configurations. The resulting joint angles and calcu¬ 
lated image pairs were used as input to the iterative algorithm given by (5.18). In ad¬ 
dition, a preliminary guess of ^ = (1.1,1.1,1.1, —.9, .1, .1, .11, .4, .11, .1, .1, .06, .06) T 
was provided and S 0 was fixed at 0.1 m. The algorithm was run until the ‘actual’ 
parameters were recovered to within eight decimal places. 


5.6 Summary 

A general framework for calibrating a manipulator and stereo system for performing 
the task of hand-eye coordination has been presented. The emphasis has been on 
autonomous calibration. The vision system was seen to be represented simply with 
the D-H convention, thus allowing a unification of manipulator and camera model 
notation. The vector based derivation of the columns of the closed-loop Jacobian en¬ 
abled an identifiability condition to be derived. As examples of the application of this 
condition the parallel axis problem, the length parameter scale problem, and the data 
‘persistency of excitation’ problem were all discussed. The iterative parameter search 
algorithm introduced in Chapter 3 was generalized and demonstrated in simulation. 


86 



Chapter 6 


Relevance of Robot Calibration to 
Human Motor Control 


The theory of closed-loop kinematic calibration is based on two premises that, though 
generally taken for granted in robotics, must be justified for human motor control: 

• There is need for a task space, and we must therefore learn transformations 
from our sensors to this 3-D internal representation - a representation distinct 
from other internal representations such as joint space. 

• There are structured internal models representing the sensorimotor transfor¬ 
mations. These models are constrained to reflect the physics of our arms and 
sensors, and do not take on an arbitrary form. 

The second premise is not independent of the first; the hypothesis that there is a 
task space requires a constrained parametric representation of the transformations in 
order to reasonably explain how the transformations are calibrated (re Section 2.4). 

Both of these premises are investigated in the following literature review. 

6.1 Internal Models 

The notion that there may exist an abstract internal model of the world may be 
traced back to Descartes [77]. For example, his contemplations of how a blind man 
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uses a stick for perception led him to the conclusion that we have available a “natural 
geometry”. Without knowing the length of the stick the blind man can triangulate 
to locate points in space. Descartes did not seem to confront the question of how 
our arms are calibrated to provide the required hand positions. Nevertheless, he did 
recognize the need for a geometric reasoning ability, or more specifically an internal 
model structure. 


6.2 Accuracy of Pointing 

6.2.1 Accuracy Without Vision 

One of the earliest systematic investigations of arm movement revealed that accuracy 
can be maintained even in the absence of vision [78]. Woodworth concluded that 
there are two phases to a movement: an impulse phase that moves the arm most of 
the way to the target, and a current controlled phase that provides final adjustments. 
During non-visually guided movement the latter phase does not exist and movement 
is executed in a feedforward manner (with respect to vision). 

To further investigate non-visually guided movements Woodworth [78] had sub¬ 
jects repeat a target movement of a certain distance. The repeat movement was 
performed (1) in a different area of the workspace, (2) with another writing imple¬ 
ment, (3) at a different speed, (4) with the wrist instead of the arm, and (5) with a 
different hand. The results show that across all of these conditions the extent of the 
movement can be reproduced to within approximately 10 percent of the extent of the 
target movement. When movements were made with different parts of the body (e.g., 
torso, foot, etc.) accuracy of repeating the extent of movement was still maintained 
to within 30 percent. Apparently the motor system is able to represent “distances,” 
independently of the sensory modality, which is consistent (but not necessarily the 
only interpretation) with the hypothesis that an internal 3-D representation of the 
world exists. 

Woodworth [78] goes on to argue that the sensation of distance is not purely 
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Figure 6-1: A subject’s attempt at drawing equal overlapping triangles with his eyes 
closed (actual size). The movement was relatively slow, ?==:500 ms per segment. 

derived from the endpoint positions in space. This is because, although the extent 
of movement can be accurately controlled, the absolute positioning is less accurately 
controlled. For example, in drawing a triangle repeatedly without vision a subject 
replicates the lengths of the sides well, but the absolute location of the vertices drift 
substantially [78]. This experiment is easy to repeat; see Figure 6-1. 

Fitts [79] also noticed that movement accuracy depends upon other variables than 
the target position: the initial position, the velocity and the extent of movement all 
influence accuracy of hitting a target. It may be speculated from these studies that 
velocity and reafferent force information are used in addition to position information. 
This conclusion is in line with present knowledge of proprioceptors. Joint capsule 
sensors seem not as important in providing position sense as spindles, and spindles 
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provide only velocity and relative length changes in the muscle [4]. Thus, the pri¬ 
mary proprioceptive sense must be integrated to provide position. Such integration is 
sensitive to initial conditions, and is prone to drift. 

6.2.2 Accuracy and Head Movement 

Fitts [79] found that blind pointing accuracy to targets initially seen visually decreased 
with the angular displacement away from straight ahead (see review [80]). Biguer 
[81] extended these results to show that the inaccuracy was principally due to the eye 
rotating laterally more than 30 degrees. If subjects are allowed to move their head 
pointing accuracy changed in the direction of the head movement. These results only 
reveal that our joint angle sensors are accurate within limited ranges - not surprising. 


6.3 Proprioceptive Distortions 

6.3.1 Normal Proprioceptive Distortions 

Other work has revealed that the position sense inaccuracies are not just due to sensor 
limitations. There are systematic distortions in our perception of distance measured 
by our arms. For example, forced choice studies reveal that distances sensed in the 
horizontal plane are perceived longer in the radial direction than across the chest 
[82, page 31-25][83]. Such global distortions cannot be accounted for by a purely 
memory-based account of the sensory transformations, and is more symptomatic of 
inaccurate structured models. Unfortunately, these results are difficult to interpret 
because distance perception is affected by duration, force and ease of movement [78]. 

6.3.2 Vibration Induced Proprioceptive Distortions 

Vibratory stimulation of the arm muscles can produce illusions of movement - presum¬ 
ably from primary spindle excitation. These illusions are so powerful that movement 
can be perceived to go beyond joint limits and distort parts of the body that the 
hand is touching. More importantly, such arm muscle vibrations can also create vi- 
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(a) 




Figure 6-2: (a) Transformations learned by strict association, (b) Transformations 
with an abstract intermediate reference frame. 

sual distortions that are most consistent with a hypothesis that a 3-D body-centered 
reference frame exists and has rotated [84]. 

6.3.3 Prism Glasses and Adaptation 

The most studied perceptual distortion is that caused by wearing prism glasses (with 
lateral displacing, rotating, or inverting lenses) [29] [85] [86] [87] [88]. In a typical ex¬ 
periment subjects wear prism glasses and are asked to adapt to manipulating objects 
with one hand under the distorted view. Complete adaption usually occurs if enough 
time is allowed (hours or days). Two aspects of the recovery are relevant. 

(1) The adapted ability generalizes from the learned set of target points. This 
global generalization supports the notion that the sensory transformations are con¬ 
strained, as opposed to being built up by strict association. It is not clear, though, 
whether the constraint is as simple as linear interpolation between neighboring points 
in space [89], or whether the transformations are due to parametric adaption of global 
models. 

(2) The other aspect of the recovery from prism distortion is that inter-manual 
transfer can occur. Inter-manual transfer refers to the ability to immediately achieve 
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accurate hand-eye coordination with the hand that was not exposed during the adap¬ 
tion period. The inter-manual transfer depends on various factors present during 
adaption: whether the eye tracks the hand or visa versa, whether the whole arm is 
seen or not, and whether neck movement is permitted or not [86] [87]. However these 
dependencies are explained, the important conclusion is that inter-manual transfer 
can occur. If inter-manual transfer never occurred then the simplest interpretation 
would be that no intermediate coordinate system is used, and direct transformations 
between the sensors are learned (Figure 6-2a). The finding that transfer does oc¬ 
cur suggests that there are separate modules for the left arm, the right arm and the 
head-eye system, and that each module transforms the sensor-space information into 
a common task space — see Figure 6-2b (there is a possibility here that task space 
is identical to joint space or retinotopic space). If adaption occurs in the head-eye 
module then inter-manual transfer is observed. If adaption erroneously occurs in the 
arm module then inter-manual transfer is not observed. 


6.3.4 Teleoperation as a Proprioceptive Distortion 

Another approach to studying adaption is to change the arm’s geometry. Though 
this cannot be done easily, an equivalent effect can be achieved by having a human 
operator control a teleoperated robot [90] [91] in a master-slave set-up (the master 
being an exo-skeleton fit on the human arm). Distortion in the kinematics could be 
systematically introduced into the slave manipulator. Variables such as whether the 
shoulder joints intersect, or the link lengths change could be controlled. As with the 
prism experiments generalization and inter-manual transfer could be studied. 

Unfortunately, the teleoperation literature does not contain such studies. There 
are reports that much training is necessary (months) [92]. A few groups have inves¬ 
tigated the problem of using differing geometry master and slave robots. A so called 
“incompatibility index” has been developed to measure the tolerated mechanical dif¬ 
ferences between arms [93]. Vertut [90] has found that up to 30 degree orientation 
mismatch is tolerated by an operator. 

I have conducted informal studies of subjects learning to control an X-Y point 
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on an oscilloscope screen with the elbow angle driving the X-movement and the 
shoulder angle driving the Y-movement. The hypothesis that the adaption is highly 
structured predicts that learning fast and accurate hand-eye (oscilloscope-eye) coor¬ 
dination in one region of the workspace should generalize immediately to the rest of 
the workspace. The two subjects tested took over three hours to learn to accurately 
move in one quarter of the workspace. By this time they could move “ballistically” 
to target points. In spite of this long adaption period, the learning immediately gen¬ 
eralized to the rest of the workspace - as predicted by the constrained representation 
hypothesis. These findings are only preliminary, and must be repeated with more 
degrees of freedom involved, and less dramatic changes in the kinematics. 


6.4 Mental Imagery 

6.4.1 Reaction Times to Recognize Rotated Objects 

The notion of an abstract 3-D reference frame is a premise of the work on mental im¬ 
agery [94]. A typical result from these studies is that the reaction time in recognizing 
a rotated object is proportional to the object’s degree of orientation away from the 
target object. The conclusion is that the object is incrementally rotated in an abstract 
3-D reference frame until it aligns with the target object. Though plausible, this abil¬ 
ity to reason spatially does not prove that a 3-D representation exists. It would be 
adequate to explain these results with an internal joint-space representation. In fact, 
in robotics joint space (or more generally, configuration space [95]) is often used for 
planning. Further, when mental rotation experiments are repeated with pictures of 
hands used as stimuli, instead of arbitrary objects, the reaction times, oddly enough, 
obey arm anatomical joint limits [96]. 

More conclusive evidence for an internal 3-D reference frame comes from studies 
in which subjects recognized letters drawn on various surfaces of their bodies. The 
letters - half mirror reversed - were drawn while the body parts were in different 
configurations. The findings turn out to be consistent with the idea that we do 
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maintain one or more internal reference frames [97]. 


6.4.2 Structure from Motion 

One of the most persuasive arguments for representing the world in a 3-D task space 
comes from our ability to recover a rigid body’s structure from the motion of moving 
points [98] [99]. Surprisingly, we can recognize a rigid object by only viewing random 
points projected from a moving object onto a screen, in spite of distance not being 
preserved between projections of points into different views. Ullman [99] argues that 
this ability requires an assumption of the object being a rigid body in the 3-D world. 
Thus, again, it seems that the motor system can represent “distance” in the 3-D 
world. 

A similar situation is encountered in joint space. The distance between points 
in joint space is not preserved as the points are rigidly moved in the world. To see 
this you only need to let one point be at your shoulder and the other at your hand. 
Movement of only the shoulder joint leaves the the distance between these two point 
fixed in space, but it does change the Euclidean distance calculated in joint space. 
Thus, none of the natural (intrinsic) sensor coordinates provide a convenient way in 
which to reason about metrical relationships in our 3-D world. 


6.5 Neurophysiological Evidence 

6.5.1 Direction Coding in the Cortex 

There are recent findings that indicate that the firing of populations of neurons in 
the cortex encode the direction of hand movement - independently of the arm con¬ 
figuration [100] or load [101]. While encouraging, these results are hard to interpret. 
The finding that cells fire independently of the actual arm configuration or load only 
indicates that there is an internal representation. It does not imply that this repre¬ 
sentation is abstract from one of the intrinsic sensor spaces (e.g., retinotopic space). 
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6.5.2 Muscles and Adaption 

It is sometimes possible to explain sensorimotor adaption by muscle or muscle spindle 
changes (e.g., see [102] and page 189 in [4]). Indeed, active movements seem to be 
necessary for adaption to prism glasses [29]. While these findings are of eventual 
interest in explaining the implementation of parametric adaption, they do not explain 
how the amounts by which to adjust the muscle parameters are calculated. Further, 
such muscle or spindle adaption can only accommodate simple changes such as joint 
angle offsets or scales (for example, changing the neck joint angle offset by adjusting 
the muscle spindles could compensate for lateral displacement prism glasses). These 
muscle mechanisms cannot compensate for errors in the knowledge of link lengths or 
joint axis orientations. 


6.6 Discussion 

In summary, the psychophysical data are not inconsistent with the hypothesis that 
we can transform sensory information into an abstract 3D task space. Further, from 
a functional stand point such a 3D task space is attractive: it enables simplicity in 
representing and recognizing objects in the world. One contribution of this thesis has 
been to show how structured models may be used to learn transformations to a 3D 
task space. Another contribution has been to show how hard this learning process 
really is, even with simplified robotics models. This difficulty of learning must be 
carefully weighed against the above simplicity arguments for using a 3D task space 
in the first place. 
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Chapter 7 


Movement Control Theories 


7.1 Introduction 

Human arm movements are surprisingly stable, fast, and accurate considering the slow 
neural hardware and large feedback propagation delays. As mentioned in Chapter 1 
this performance is partly due to pre-planned feedforward motor drive. We regularly 
make ballistic movements without visual guidance, and animals are found to still 
execute normal movement after proprioceptive deafferentation [103]. In the previous 
chapters the focus was on the learning of feedforward kinematic information. In the 
following chapters dynamic feedforward compensation is studied. The issue is not 
how a feedforward dynamic model is learned, but whether such a model exists, and 
what form it takes. 

Another reason for the high performance of the human arm is muscle dynamics. 
Muscle has a very large dynamic force range. It can produce forces per unit weight 
exceeding current robotics actuator technology, while it can still produce fine enough 
forces to control the hand of an eye surgeon. Muscle can also be extremely fast; for 
example, contraction frequencies of 120 Hz have been observed in some insects [104]. 
But most importantly, muscle has built in compliance. It has a short range high 
stiffness, which may function in impact situations such as running [105], and a longer 
range low stiffness [105], which may allow compliant motion control without excessive 
active feedback. 
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A series of experiments are described with the goal of characterizing the control 
system used in arm movement. Feedforward compensation and muscle properties will 
be shown to be of particular importance. 


7.2 Non-Linear Operator Description of Arm Con¬ 
trol 

The arm control system may take on the general form depicted in Figure 7-1, where 
it is assumed that the task is to track a joint angle trajectory. The input co mman d is 
the desired joint angle trajectory 9d(t) and the actual output trajectory is 9{t). The 
non-linear operator A represents the feedforward dynamics. It computes the feedfor¬ 
ward torque Td(t) = A9d(t). (The word operator is used here in the mathematical 
sense defined in [106] [10]. An operator maps one whole time sequence onto another 
whole time sequence. In contrast to a function, an operator may represent systems 
with memory.) The non-linear operator 71 represents the feedback system driven by 
muscle spindles (or other proprioception) with a sensitivity controlled by the gamma 
drive. Finally, the non-linear operator V represents the musculo-skeletal plant being 
controlled. Though the following treatment will be general, it is helpful for under¬ 
standing to keep in mind the example of a proportional-plus-derivative controller 
acting on a second order linear system. In this example V~ x = Id 2 /dt 2 + Bd/dt + K, 
H. = K v d/dt + K p and >1=0. 

In the absence of external perturbations r, inspection of Figure 7-1 yields the 
general input-output relation of: 

(1 + VH)9 = (VA + VK)9 d (7.1) 

If A inverts the musculo-skeletal plant dynamics (i.e., VA = 1) then perfect feedfor¬ 
ward control is achieved: 9 = 6 d (except for input signals zeroed by 1 -f VTV). Also, 
if the open-loop operator W has a much larger effect than either VA or the unity 
operator 1 on any signal (i.e., the ‘gain’ [10][106] of W, denoted \V7l\, is large), 
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Figure 7-1: General controller for a joint angle trajectory 0(£). 

then high gain feedback command following is achieved: $ « 9j, (except for input 
signals zeroed by VTV). This last item should be used with consideration of the small 
gains theorem [10] [106]: if the open-loop gain \V7l\ is less than unity, then stability 
is always guaranteed. See Section 1.1.1. 

7.3 Closed-Loop Dynamics and a Perturbation 
Model 

With human subjects, aside from measuring EMG, the motor system may only be 
studied by injecting external torques r(£) and measuring the joint angles 9{t). None 
of the other quantities in Figure 7-1 can be measured. In this section a model of 
relationship between total torque inputs and joint angle displacement will be devel¬ 
oped for the motor system governing the elbow joint. Then a perturbation model 
will be developed to relate externally applied perturbations to recorded joint angle 
movements. 


99 



7.3.1 Model of Closed-Loop Single Joint Dynamics 

The operator expressing the closed-loop dynamics relating total torque inputs to joint 
angle 9{t) is easily derived to be V c i such that: 


(i + vnypa = v ( 7 . 2 ) 

We now derive a more detailed form of V c i. First transform Figure 7-1 into the 
functionally equivalent form of Figure 7-2: 



Figure 7-2: Transformed equivalent diagram to Figure 7-1. 

where A! = A.+ TZ, and r„(t) is considered to be the commanded (voluntary) torque 
input to the feedback controller. Now, include a muscle model as follows. 

7.3.2 Model of Muscle 

Assuming negligible serial tendon compliance, the single joint musculo-skeletal dy¬ 
namics 'P may be represented by the following general circuit diagram: 
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Figure 7-3: Musculo-skeletal dynamics V. 

where I is the inertia of the forearm, e is the muscle torque generation mechanism 
and z(9(t), 0(t), t ) represents the intrinsic parallel force dynamics of the muscles. Of 
course, multiple parallel muscles may act on the arm, but their effect can always be 
reduced to a circuit equivalent to Figure 7-3. Symbolically, Figure 7-3 may be written 
as: 

e = I9(t) + z(9(t),9(t),t) (7.3) 

Defining the new operators Z and X such that: 

Z$(t) = z(9(t),9(t),t) (7.4) 

I9(t) = I9(t) (7.5) 

the muscle dynamics may be explicitly included in Figure 7-2 to give Figure 7-4: 
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Figure 7-4: Closed-loop dynamics with muscle model explicitly shown. 

7.3.3 Equivalent Feedback Model of the Joint Mechanical 
Properties 

A further transformation of the diagram in Figure 7-4 gives the functionally equivalent 
Figure 7-5: 



Figure 7-5: Equivalent feedback model of the joint mechanical properties, 
or symbolically Figure 7-5 is: 

0 = I0(t) + (Z + K)0(t)-T 9 (t)-T(t) (7.6) 
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This last loop transformation provides an important insight that will be used repeat¬ 
edly. That is, the intrinsic muscle dynamics can be equivalently built with an active 
feedback system, and effectively provide feedback control. This is the idea behind the 
equilibrium control theories of [28], explaining experiments in which active feedback 
was removed. We will refer to the effective feedback system and will mean that of 
Figure 7-5 — generated by both intrinsic muscle properties Z and active feedback 71. 

Just as the mechanical compliance of the muscles may be built from an equivalent 
active feedback controller, the active feedback can sometimes be realized by a phys¬ 
ically equivalent passive device (e.g., Figure 7-5 can be realized with an equivalent 
physical circuit of the form shown in Figure 7-3 with / = z+r substituted for z, where 
r is defined from the operator 7 Z as described below). Colgate [11] elaborates on this 
with regard to stability. We will also use this physical equivalence notion implicitly 
when we talk about mechanical compliance. The term compliance is typically applied 
to passive mechanical devices, but as an actively controlled feedback system can often 
be built with an equivalent passive device we can say that this active system has an 
effective compliance , or compliance for short. 

7.3.4 Model of Reflexes 

To proceed, assume that the proprioceptive feedback only provides velocity and po¬ 
sition information. Further, let us neglect the effect of feedback delay (this will 
be discussed further in Section 7.6). Then, in general the non-linear feedback is 
719 = r(0(t),0(t),t), and (7.6) becomes: 

0 = I$(t) + r(6(t), 9(t), t ) + z(0(f), *(*), t) - r v (t) - r(t) (7.7) 

7.3.5 External Torques 

Considering the perturbation torque r(f) to be the input, and adding the possibility 
of gravity (7.7) becomes: 
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(7.8) 


r(f) = IO(t) + f(0(t),0(t),t) — r v {t) — mgc sin 6{t) 

where 7, m and c are the forearm inertia, mass and center of mass respectively, 

— r(9(t),0(t),t) + z(0(t),0(t),t), and g = 9.8 m/s 2 . Notice that 
f — t v is the net muscle torque acting on the arm inertia (including parallel connective 
tissue forces, active feedback, and intrinsic mechanical properties). The gravity term 
was derived by assuming that the arm is moving in the vertical plane, and 6 is defined 
to be zero at vertical as described in Section 9.1. 

7.3.6 Perturbation Model Used in Experiments 

We are now ready to address the question of how to experimentally study the closed 
loop dynamics V c i represented in (7.8). As the internal forces (e.g., r v (t)) that act 
in addition to the known external forces cannot be controlled or measured, only the 
local behavior of V e i about a nominal trajectory 9 0 (t) may be studied. Applying an 
external perturbation r(<) during execution of the nominal trajectory gives a new tra¬ 
jectory that is related to the nominal trajectory through the dynamic behavior of the 
arm. The arm dynamics can be studied by estimating the relationship between the 
applied perturbations and the deviations from the nominal trajectory. Specifically, 
if the applied perturbations are small enough, the closed-loop dynamics may be ap¬ 
proximated by the differential behavior about the operating point (0 o (t),T O (f),r„(<)), 
which is: 


Ar(f) = I(t)A0(t) + B(t)A6(t) + K(t)A0(t) (7.9) 

A 9{t) = 9(t)-0 o (t) 

A T (t) = r(t) - r 0 (f) 
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where the joint inertia I(t ) = I + df/dO , the joint viscosity B(t) = df /d9 and the 
joint stiffness K(t) = df/dO — mg c cos 0 o (t). The second term of I(t) should be zero 
as it could only result from high bandwidth force feedback. Also, the negative gravity 
contribution to the stiffness turns out to be relatively small (see Figure 10-16 in the 
results, Chapter 9). 0 o and r 0 are the angles and torques measured during an average 
unperturbed movement. r 0 may be non-zero. 

In the sense that the parameters /, B , and K relate force input to velocity output, 
they may be referred to as mechanical impedance parameters. The relation between 
force input and position output may be referred to as complex stiffness , with compli¬ 
ance being the inverse of this relation. 


7.4 Control Hypotheses and Predictions 

In subsequent chapters experiments to estimate the perturbation model (7.9) are 
described. In these experiments subjects perform the simple task of moving their arm 
between two targets, and perturbations are applied to estimate the dynamics. For the 
present, assume that the locally Unear behavior of "P c / may be estimated from these 
perturbations, and that it may be adequately modeled by the time-varying second 
order system (7.9). Now, consider what this experimentaUy determined perturbation 
model may reveal about the modules A, H and V in Figure 7-1. 

7.4.1 Feedforward Versus Feedback Control 

First, focus on the role of A, the feedforward control. One extreme hypothesis holds 
that A compensates for the plant dynamics entirely (i.e., A = V -1 , so 0 = do). 
This feedforward compensation is non-trivial (i.e., not a Unear system) for multi- 
link arms [21] or even for a single joint operating against gravity. In spite of the 
difficulty of learning and computing 'P~ x , this hypothesis of feedforward compensation 
is consistent with experimental evidence. SpecificaUy, trajectories are found to be 
invariant across load, gravity and speed conditions (caU this trajectory invariance ). 
That is, for a given movement the trajectories are the same in spite of differing inertial 
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loads or gravity [107] [108], and even take on the same form for different speeds once 
normalized by movement time [107]. 

A competing explanation for trajectory invariance is that the closed-loop system 
V c i is made to operate as an effective high gain feedback controller (| V7l\ » 1, so 
9 = Qq). Non-linear link interactions [21] or gravity are treated as disturbances to be 
servoed out by the high gain feedback. As shown in Figure 7-5, this may be imple¬ 
mented either by active feedback, 7 Z, or by altering the intrinsic muscle properties, 
Z. Both of these possibilities are unlikely. With regard to reflexes, Bizzi et al. [103] 
provide evidence showing that reflex gains are low. This is not surprising, as delays 
of more than 25 ms would make high gain feedback unstable. With regard to muscle 
stiffness, for normal speed movements the co-contraction necessary to produce high 
intrinsic stiffness is rarely observed [109][110] — tri-phasic, or bi-phasic alternating 
muscle bursts are usually observed. Thus, neither of these implementations are ac¬ 
ceptable, and an account of trajectory invariance requires some form of feedforward 
control. 

7.4.2 Compliance Control Hypothesis 

A more likely scenario is that feedforward control is used, but it is not always ac¬ 
curate enough, so low gain feedback is used to assist when necessary. In this way 
trajectory invariance may be accounted for, and feedback 71 and intrinsic muscle 
properties Z may still have a role. More generally that role would be to control 
the mechanical compliance of the arm. Specific examples of compliance control have 
been suggested previously. For example, the reflex action of feedback may make the 
system dynamics look more like a linear compliance, thus simplifying the job of the 
feedforward controller [111]. Another suggestion is that low levels of muscle stiffness 
may provide stability against unexpected disturbances [28]. Finally, as mentioned 
in Chapter 1, the control goal of having a high mechanical compliance is advanta¬ 
geous for constrained motion, particularly when the constraint is unpredictable or 
unexpected (e.g., to avoid jamming, excessive contact force, damage or injury). 

Thus, for the movements studied in this thesis a reasonable control hypothesis to 
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test is that the arm’s mechanical compliance is maximized, while still maintaining 
the specified task accuracy. That is, the arm’s effective stiffness (i.e., restoring forces 
to trajectory errors — with the motor system viewed as the feedback controller in 
Figure 7-5) should be low, but high enough to assure that trajectory errors larger 
than the specified task accuracy are removed. As imperfect feedforward control or 
external perturbations will inevitably cause trajectories errors, the effective stiffness 
must be positive at all times. In summary, the first prediction for the perturbation 
experiments is: 

(Pi) In so far as the arm’s stiffness can be adjusted for a particular task, the stiffness 
should be as low as the task accuracy allows. In the task of pointing between targets 
accuracy during movement is not demanded. Thus, the stiffness should be low during 
movement (to provide compliance), and higher at the targets (to ensure accuracy). 

The mechanism for controlling arm compliance may rely on both intrinsic muscle 
properties and reflex activity. For example, the intrinsic velocity dependent dynamics 
of muscle could make lowering of the stiffness during movement automatic (this is 
discussed further in Sections 7.5 and 11.3). In addition, reflexes are known to ad¬ 
just dramatically during particular tasks. For example, in walking there is a strong 
anterior tibial H-reflex during stance, but little reflex response in the swing phase 
[ 112 ]. 

7.4.3 Limited (Static) Feedforward Control Hypothesis 

The computational complexity of feedforward control naturally leads to the question 
of how much feedforward control is really necessary. Consider the hypothesis that only 
the instantaneous static component of the dynamics V c i is inverted by the feedforward 
A! in Figure 7-2. That is, the feedforward torque t v is computed by (7.8) with 
9 = 0 = t = 0: 


T v{t) = f(0 d (t),0,t) -mgcsm0 d (t) (7-10) 
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Under this hypothesis, the desired trajectory 9 d (t ) can be viewed as an equilibrium 
trajectory [28], in the sense that it is only achieved when the arm comes to rest. 
During movement of the arm there is always a trajectory error, due to the non-static 
components of V c i. Note that this is feedforward control in the sense that it requires 
knowledge of the static properties of muscles. The behavior of the actual trajectory 
9(t) driven by the static feedforward command (7.10) is governed by: 


0 = I9(t) + f{9(t),9(t),t)-f(9 d (t),0,t) 

—mgc( sin 9{t) — sin 9 d (t)) (7.11) 

Now consider the consequences of this static feedforward hypothesis for attempting 
to generating two identical but different speed movements. That is, suppose (7.11) 
describes a particular movement 9(t) specified by the equilibrium trajectory 9 d (t), 
and suppose a faster speed movement 9f(t) is attempted under this hypothesis — 
that is by simply time scaling the equilibrium profile by a factor of r. Thus, the new 
desired trajectory is 9 d (rt). Substituting this into (7.11), and replacing t by t/r gives 
the equation governing the behavior of the faster movement #/(<): 


0 = r 2 I9 f (t/r) +f(9f{t/r),r9 f (t/r),t/r)~ f(9 d (t),0,t/r) 

—mgc (sin 9f(t/r) - sin 9 d {tj) (7.12) 

where the time scaled faster movement 9f(t/r ) should be close to 9{t). 

To see how different 9f(t/r) and 9{t) actually are make the following simplifica¬ 
tions: (1) assume the dynamics / are time invariant (or at least they change at a rate 
r faster for the faster movement), and (2) assume that the dynamics are the same 
for both movement speeds (later we will see how scaling the dynamics may change 
things). Now subtract (7.11) from (7.12) to give: 
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(7.13) 


0 = I(r 2 i),(t/r) - «'(()) + /(<!,((/r), r«,(t/r)) - /(0(f), 0(f)) 

— mgc(sm0f(t/r) — sin0(t)) 

When the movements just start (without loss of generality, let that be at 0 = 0) only 
the acceleration terms are significant. Thus, (7.13) implies r 2 0f(t/r) = 9{t). That 
is, there is a factor of r 2 less acceleration in the time scaled faster movement at the 
beginning the movement. Equivalently, the scaled faster movement looks like it has 
r 2 times more inertia in its dynamics. Once significant velocity is reached predictions 
are harder to make without simulations, but the trajectories do diverge increasingly. 

Alternatively, the dynamics might be adjusted to compensate for an increase in 
movement speed. Consider the case where the dynamics are linear and of second 
order. In this case, increasing the stiffness by a factor of r 2 and the damping by a 
factor of r is the appropriate adjustment to have the faster movement 0f(t) have a 
time scaled trajectory 9f(t/r) equal to the slower trajectory 0(t) [107][113] (assuming 
no gravity is present). To prove this, let f{9f{r),9 f (t),t) = B9 f (t ) + K0 f (t), and 
g = 0 in (7.12), and divide through by r 2 to give: 


0 = I9f(tlr) + r- 1 B9 f {tlr) + r- 2 K{0 f {tlr)-9 d {t)) (7.14) 

Letting B = tB and K = r 2 K eliminates r, and thus makes the time scaled solution 
9f(t/ r ) identical for all speed movements. 

In summary, with regard to different speed movements the static feedforward 
hypothesis implies the following prediction for the perturbation experiments: 

(P2) Either the time scaled trajectories must differ for different speed movements, or 
the dynamics must be scaled. If the dynamics are approximated by a second order 
linear system, this dynamic scaling requires that the stiffness increases by a factor of 
r 2 and damping increases by a factor ofr, where r is the movement speed increase. 
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Next consider the influence of gravity on a control system operated under the 
static feedforward hypothesis. The prediction for the experiments is simply: 

(PS) The absence or presence of gravity should not affect the stiffness, damping or 
trajectory. It should be compensated for by the hypothesized static feedforward system. 
Of course, if the actual trajectory differs significantly from the desired trajectory the 
gravity compensation might be inappropriate. See (7.11). This deviation of the actual 
trajectory from the desired trajectory should be small for slow speed movements, but 
increase with increased speed of movement — as discussed above. 

It is possible that gravity compensation is achieved by feedback instead of feed¬ 
forward control. For example, there could be an integrator in the feedback controller 
7£, eliminating the steady-state positioning error (e.g., PID control). In this case, the 
linearized dynamics would have to be at least third order. As second order dynam¬ 
ics adequately account for the variance in previous posture perturbation experiments 
(see review in Chapter 8) this is not expected to be a possibility. 

Finally, with respect to non-static changes in the dynamics the static feedforward 
hypothesis does not predict compensation. For example, 

(P4) When the dynamics are changed by an external viscous damping load, the static 
feedforward compensation should not adjust for this change, and the trajectories should 
thus differ from the unloaded trajectories. 

Earlier studies of whole arm movements indicate that damping loads do change the 
trajectories. Using the instrumentation techniques and methods of [108][107] we stud¬ 
ied whole arm sagittal plane movements with a damping load applied to the elbow 
joint. Figure 7-6 shows typical normal speed trajectories made with and without 
the damping load of 1.5 Nm/rad/s. The figures show three typical movements be¬ 
tween targets after 20 practice trials. Movements recorded immediately after changing 
the damping are very similar, showing little adaption. Particularly striking is the in¬ 
creased difference between the upward (dotted lines) and downward movements (solid 
lines). See also [114]. 

The finding that a damping load does affect the arm’s trajectories is to be con¬ 
trasted with the finding that whole arm trajectories are unaffected by a mass load 
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Figure 7-6: Whole arm joint-space trajectories made (a) without and (b) with a 
damping load of 1.5 Nm/rad/s acting about the elbow joint. Subjects were given time 
to adapt to the load changes (though little adaption occurred), and in the unloaded 
case the subjects still carried the mass of the damper. Three upward (solid) and three 
downwards movements (dotted) are shown. Each movement took approximately 500 
ms. 


carried in the hand [107]. If the results in [107] are correct, there must at least be 
non-static compensation for inertial loads. 

Approximate estimate of desired trajectory 

Under the rough approximation that the muscle dynamics are linear it is possi¬ 
ble to back calculate the desired equilibrium trajectory, given a static feedforward 
hypothesis: 


9desired = 9 + l/K(t)[I(t)9 -f- B(t)0 — mgc sin 0] (7-15) 

The difference (0 — 9desired) provides an indication of the positioning accuracy of the 
feedforward compensation. Although the dynamics may be non-linear (though see 
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[Ill] for the possible linearizing effect of reflexes), this calculation is approximately 
correct if the 9, 9 and (9 — 9 d e tired) are respectively of similar amplitude and frequency 
content to the A 9, A 9 and A 9 perturbations used to estimate J, B , and K. 

7.4.4 Dynamic Scaling Feedforward Control Hypothesis 

Another possible method to simplify feedforward control is to learn accurate feed¬ 
forward commands for only one slow speed movement and scale up the torques and 
dynamics to achieve the same movement at higher speeds [107] [113]. Assuming mus¬ 
cle is adequately approximated by a linear system it is possible to predict under this 
dynamic scaling hypothesis that: 

(P5) The stiffness should scale up by a factor ofr 2 , and the viscosity should scale up 
by a factor of r, where r is the speed increase factor. 

The weakness with this argument is, again, that the dynamics may not be linear. 

For single joint movements executed under this dynamic scaling hypothesis the 
required stiffness and viscosity scaling keeps the damping ratio fixed. In recent pos¬ 
ture experiments we have found that the damping ratio is constant across three levels 
of voluntary muscle co-contraction [115] [116]. Scaling of the dynamics during move¬ 
ments has not yet been tested. 

7.5 Role of Muscle Dynamics 

The measured perturbation dynamics are a function of active feedback, inertia, con¬ 
nective tissue, and muscle. To assess the possible contribution of muscle the following 
review is provided. 

7.5.1 Frequency Range of Interest 

First, consider what aspects of muscle dynamics are not relevant. For the uncon¬ 
strained movements studied in this thesis, it is possible to show that above 5 Hz the 
elbow joint dynamics are dominated by the arm inertia (see Figure 10-8). Thus, only 
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muscle/tendon dynamics between 0 and 5 Hz are relevant to understanding these 
arm movements. 

7.5.2 Molecular Level Description 

On a molecular scale the mechanism of muscle force production is described well by a 
model of crossbridge bonding between sliding filaments [117]. A cycle of crossbridge 
connection and disconnection happens continuously. The instantaneous number of 
connected crossbridges is related to the stimulation frequency, and it is roughly pro¬ 
portional to the net muscle force. 

When a tetanized muscle fiber is subjected to an externally imposed step change 
in length there results a large transient increase in force (typically, lasting only a few 
milliseconds), and then a gradual change in force to a steady-state value [118]. The 
latter force dynamics are most relevant to the present arm movement study, (i.e., they 
are in the 0-5 Hz frequency range, above which inertia dominates). To a first ap¬ 
proximation, the transient force can be attributed to the stiffness of the crossbridges 
before they break, and the dynamics of the slower force production process can be 
attributed to the crossbridges cycling and reforming. A more detailed analysis re¬ 
quires consideration of the distributed nature of the sarcomeres within the fiber. The 
sarcomere lengths are known to be non-uniform and change even when the muscle is 
in isometric contraction (see review in [119]). 

7.5.3 Steady-State Length-Tension Curves 

The steady state force value is a function of only the final length, giving the so called 
length-tension relation. Figure 7-7 shows a length-tension plot made by Gordon et 
al. [120] from an experiment where they measured the steady state force reached in 
response to an isotonic change to a specified fixed length (some points extrapolated). 
The force change to steady state takes 50-100 ms (as in [118]), which implies that the 
slope of the length-tension curve can only be interpreted as a ‘stiffness’ when dealing 
with phenomena slower than 50-100 ms. 
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Figure 7-7: Static length-tension curve for changes of sarcomere length in a muscle 
fiber. Numbers correspond to interpretation in terms of the sliding filament theory. 
Reprinted from [120]. 



7.5.4 Isotonic Shortening 

During muscle movement the length-tension relation has no meaning. In fact, during 
the isotonic shortening of the muscle fiber in [120] the velocity is relatively constant, 
indicating that the stretching fiber presents a viscous-like load, with zero stiffness. 


7.5.5 Whole Muscle 

On the macroscopic scale similar length-tension relations are found [121] [122] [123] 
(see Figure 7-8 from [122]), although the mechanisms are more complex. Mammalian 
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Fig. 3. Tension during extension through different parts of the length range. 
Tension records from a number of different extensions similar to those 
shown in Fig. 2 have been traced (continuous line) to show how the tension 
during lengthening differs from the isometric tension in various parts of the 
length range. Positions of the ankle joint are shown below corresponding 
parts of the abscissa. 

In each case the muscle was lengthened through 6 mm at 7-2 mm/sec. 
The interrupted lines are isometric length-tension plots. 


Figure 7-8: Static length-tension curves for responses to changes in whole muscle 
length, with different rates of asynchronous stimulation. Reprinted from [122] 

muscle fibers are arranged in parallel, and attached to ligamentous tissue in a com¬ 
plex pennate structure. Changes in pennation angle can allow muscle length changes 
without significant muscle fiber length changes [124]. Also, muscle fiber recruitment 
can alter the dynamics. Muscle is innervated by many intermixed motor-neurons, 
which are commonly thought to be recruited by size [125], but can also be recruited 
in complex synergistic patterns [126]. Net muscle force production is regulated by 
recruitment and stimulation rate. At low forces recruitment is found to be the domi¬ 
nant control mechanism, and thus affects the muscle stiffness most dramatically [127]. 
Presumably, the force and dynamic stiffness should increase with an increase in the 
number of parallel muscle fibers recruited. 
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7.6 Role of Reflexes 


The experiments described in this thesis do not provide a clear distinction between 
reflex induced stiflness and intrinsic muscle stiffness. An ideal way to investigate 
the role of feedback would be to reversibly block all relevant afferent information 
(e.g., with a fusimotor block [128] [129]) during selected movements and observe the 
difference in the estimated perturbation parameters. As this is a difficult experiment 
(even in animals) we can only distinguish the feedback and intrinsic muscle properties 
from the closed-loop dynamics, V c i. This distinction is perhaps possible because of 
delay. 

If delayed feedback 71 is significant, then the closed-loop dynamics V c i should 
be a function of not only the current state (velocity and position), but also delayed 
velocity and position. Thus, to a first approximation the linearized dynamics should 
be a second order model with delay terms: 

A r(t) = I(t)A0(t) + B(t)A9(t) + K(t)A9(t) (7.16) 

+B d (t)A9(t + T) + K d (t)A9(t + T) 

where, for example, T = 25 ms. The parameters of (7.17) can be estimated, and the 
relative contribution of the delayed terms can then be tested. This will be left for 
future work. 
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Chapter 8 


Literature Review 


Many perturbation studies have been performed to infer intact arm dynamics. This 
chapter will review these studies to compare instrumentation, modeling techniques, 
and perturbation parameter estimates. With the exception of [130] the research has 
been restricted to measuring the mechanical properties during posture. 

The review is organized in terms of the instrumentation used for applying the 
perturbations and the perturbation type. This organization is perhaps the clearest as 
the particular instrumentation limits the range of perturbations possible (and in some 
cases changes the arm’s dynamics), which in turn affects the estimated perturbation 
dynamics. If the motor system was linear all perturbation methods should give the 
same results. Unfortunately, it is non-linear, and impedance estimates are known to 
depend on the perturbation amplitude, type, and frequency content. 


8.1 Passive Mechanical Perturbations 

8.1.1 Sinusoidal Position Perturbations 

A simple method of perturbing the arm is to attach it to a shaft mounted eccentrically 
on a rotating fly wheel (a rotating cam setup). In [12] (see also Zahalak [131]) 
such an apparatus was used to apply sinusoidal position perturbations of up to 22 
Hz. Perturbation amplitudes of 0.46 mm, 0.96 mm, and 4.6 mm were tested while 
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the subject maintained large mean forces (35-130 N ). The frequency dependent 
damping was inferred from the out-of-phase force responses. Damping increased 
roughly proportionally to frequency, as is expected of a second order system, with 
the exception of the behavior between 8-12 Hz. In the 8-12 Hz frequency range the 
damping surprisingly dropped below zero for the smallest amplitude perturbations, 
predicting a form of limit cycle instability. 

The main problem with this approach is that the low frequency sinusoidal stimuli 
may be tracked voluntarily. At frequencies below 6 Hz the authors found that the 
results were very variable and could not be used. As the in-phase system response 
was dominated by the arm inertia at or above 6 Hz (i.e., the in-phase force increased 
proportionally to the square of the frequency) stiffness estimates were difficult to 
make. 


8.1.2 Viscoelastic Perturbations 

Another method to assess the mechanical response of the arm is to attach it to a 
pair of springs [13] [132]. The arm is set into oscillation by off-setting the spring from 
the equilibrium. In [13] this method was used to drive the arm into a sustained 8-12 
Hz, 3 mm tremor (their predicted range of instability [12]). In [132] this method was 
used to show that the arm’s natural frequency (considered as a second order system) 
increases with the subject’s voluntary resistance to the perturbations. They also used 
a pseudorandom binary input (see below) to obtain a more detailed description of 
how the dynamics change with voluntary command change. 

8.1.3 Inertial Perturbations 

Bizzi et al. [133] used inertial loads to assess the difference in monkey head posi¬ 
tioning with and without proprioceptive feedback. They found that even without 
proprioceptive feedback a monkey could still move its inertially loaded head to a de¬ 
sired position. As a follow up to this experiment Bizzi et al. [134][135] [136] performed 
a similar experiment on monkey arms. Here they perturbed the deafferented arm 
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with a position pulse during movement [136], and found that the arm qualitatively 
returned to the original path after the pulse was over. Thus, they speculated that 
a virtual trajectory was specified as the input command, and the intrinsic muscle 
stiffness maintained the arm on that trajectory. No stiffness estimates were made 
during the movements. 


8.2 Electromagnetic Actuators 

8.2.1 Sinusoidal Positions Perturbations 

Agarwal and Gottlieb [137] studied the position responses of the ankle joint to sinu¬ 
soidal torque inputs applied by a DC torque motor. They also found that sinusoidal 
inputs were problematic. Such inputs produced a different compliance than measured 
with random inputs (see below). They suggest that the motor system’s dynamics 
changes to adapt to each separate frequency input. 

8.2.2 Sinusoidal Force Perturbations During Movement 

Lanman [130] studied the response to sinusoidal force perturbation applied during 
movement. Only high frequency force inputs between 15 and 30 Hz were analyzed. 
He used a similar analysis technique to [12], which relies on a preliminary estimate 
of the arm inertia to estimate the total joint torque. This total joint torque was then 
compared with joint acceleration. He assumed a frequency dependent second order 
model, defining the in-phase component of the joint torque as the stiffness, and the 
out-of-phase component as the damping. The main finding was that the joint stiffness 
is proportional to (1) the load force and (2) the inverse of the velocity. For unloaded 
movements the stiffness drops during the movement. For movements loaded with a 
viscous damper the force dependency of stiffness cancels out the velocity dependency, 
and can even increase the stiffness over the resting stiffness. Lanman attributed the 
measured stiffness to short range stiffness of crossbridges (see discussion in Section 
11.3). 
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8.2.3 Step Force Perturbations 

Crossman and Goodeve [8] studied wrist rotation responses to small amplitude step 
changes in torque. Although they do not give details of the size of the perturbation 
or the nature of force actuation, they claim that a critically damped second order 
system models the responses well. 

8.2.4 Step Position Perturbations 

Mussa-Ivaldi et al. [138] used a two link direct drive manipulandum to apply pertur¬ 
bations to the whole arm supported in the horizontal plane at posture. They used a 
‘do not intervene’ instruction to the subject, and defined the stiffness as the ratio of 
force to position change after about 200 ms. They found that the two dimensional 
stiffness fields always had a characteristic maximum value oriented radially through 
the shoulder. 

Sinkjaer et al. [139] studied human ankle torque responses to step changes in 
angle. A servo controlled geared DC motor was used to provide 2° step angle changes 
with a 40 ms rise time (other amplitudes between 1° to 7° were also used). The 
‘steady state’ (after 450 ms) torque changes were used as a measure of stiffness, with 
the subjects being asked to ‘not intervene’ during the perturbations. In spite of the 
possibility of voluntary intervention, they found that the stiffnesses measured at low 
background torques were comparable with the stiffnesses measured for the ankle joint 
with pseudorandom inputs [140] [141]. However, at higher background torques they 
found almost constant stiffness, contrary to [141]. 

An interesting innovation in [139] is that they were able to measure active muscle 
stiffness without reflexes. The authors found that by stimulating the deep peroneal 
nerve transcutaneously they could activate the anterior tibial muscle to both maintain 
a fixed mean torque and abolish reflex or voluntary input. These intrinsic muscle 
stiffnesses were found to be only slightly lower than the stiffnesses measured with the 
reflex intact. 
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8.2.5 Gaussian Random Torque Perturbations 

Agarwal and Gottlieb [140] studied the position responses of the ankle joint to random 
Gaussian noise torque input. Band limited (2 to 30 Hz useful power) Gaussian white 
noise was applied with a DC torque motor (the authors neglected to mention the 
amplitude (i.e., standard deviation) of the input). They found that, in contrast 
to their earlier sinusoidal perturbation studies, the response could be adequately 
modeled by a linear second order system (with coherence of about 90%). 

The ankle joint stiffness was found to increase with mean voluntary torque pro¬ 
duction (ranging from 17 Nm/rad at zero torque to 45 Nm at 2 Nm torque) and was 
angle dependent over a 24° range (ranging from 15 Nm to 30 Nm for the zero vol¬ 
untary torque level). Using the estimated inertia value of approximately 0.02 kgm 2 
(which includes the inertia of the apparatus) these stiffness estimates give natural 
frequencies between 4 and 7 Hz. 

Yosef and Inbar [142] also used Gaussian noise torque inputs to measure joint 
impedance. They used a torque motor to apply input forces to the elbow joint, but 
complicated matters by not measuring the motor torque. Thus they had to estimate 
the torque from the motor dynamics. They found that the estimated joint inertia can 
change and attributed these changes to Golgi organ tendon force feedback. Another 
explanation is that their estimated torques were not accurate. 

8.2.6 Pseudorandom Binary Torque Perturbations 

Soechting et al. [143] (see also [144] [132] [145] for related work by these authors) used 
pseudorandom binary sequence (PRBS) torque inputs to characterize the EMG re¬ 
sponse of the elbow joint during movement. They used an amplitude of 8 Nm, and 20 
ms elements in the PRBS. A novel time-varying impulse response identification tech¬ 
nique was used to characterize the force to EMG response. Although they recorded 
angular position information, they unfortunately did not report any force-position 
transfer functions. 
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8.3 Hydraulic Actuators 


Hunter and Kearney [146] [147] used angular pseudorandom binary sequence (PRBS) 
inputs to characterize the human ankle joint torque response. They drove a position 
servoed electro-hydraulic actuator with a 200 Hz PRBS command to produce an 
input power spectrum flat from 0 to 25 Hz. They found that a second order model 
adequately accounted for the response, but that the model parameters depended upon 
the perturbation amplitude. 

Weiss, Hunter and Kearney [148] [149] [141] extended the results of [140]. They 
found that there is a position dependence of the joint stiffness over the entire range 
of motion, and that there is a linear relation between mean voluntary torque and 
stiffness for each position. 

Damping ratio’s of 0.3 to 0.4 were measured (the same as for the random inputs 
in [140]), which should be contrasted with the damping ratio’s of 0.8 to 1.0 measured 
in pulse perturbation experiments (e.g., [8]). See Section 11.1 for further remarks 
concerning this discrepancy. 

They also found that at high voluntary torque levels the estimated inertia in¬ 
creased. This anomaly is likely due to the muscle stiffness increasing so high that the 
actuator could not produce sufficient power at frequencies where the inertia domi¬ 
nated. Thus, for high force levels frequencies of up to 100 Hz are required for proper 
system identification. 

8.4 Pneumatic Actuators 

8.4.1 Pseudorandom Binary Torque Perturbations 

We have reported preliminary results with a wrist-mounted pneumatic thruster (air- 
jet) [116][150][151]. This one dimensional 400 g airjet applies ±4 N binary forces 
to the wrist, at frequencies up to 75 Hz. The airjet switches the flow direction by 
a fluidic property called the Coanda effect. The current limitations of this airjet 
are that it can only produce binary force sequences, and that the magnitude of the 
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force cannot exceed 4 N. The advantages are that it is unencumbering, it changes 
the arm dynamics minimally, and it has a high bandwidth. We intend to develop a 
3-dimensional version similar to Colgate’s airjet described below [152]. 

8.4.2 Graded 3-Dimensional Torque Perturbations 

Earlier wrist mounted airjets were designed by Colgate [152] and Murray [153]. Both 
of these designs incorporated the additional constraint that the airjet had to produce 
graded forces. The airjet in [153] produces 8 N pulses with a 40 ms rise time. It is 
of low bandwidth, but has three states: +8 N , -8 N , and 0 N. The airjet device in 
[152] produces controllable 3-dimensional forces at up to 8 Hz. The airjet’s open-loop 
performance (i.e., just producing binary forces) could exceed 25 Hz. Both of these 
designs were based on a spool-valve switching mechanism. 
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Chapter 9 


Time-Varying Identification 
Method 


The objective of the experiments is to find the time-varying parameters /(t), B(t) 
and K(t) in (7.9) while the subject is executing a fixed movement pattern, specified 
by B 0 (t). It should be emphasized that these impedance parameters are only valid 
for the operating conditions under which they are measured. In no way is the system 
assumed to be linear. This perturbation model describes the locally linear behavior 
of the closed-loop dynamics V c i. 


9.1 The Task 

In order to estimate time-varying dynamics many repeated movements are required to 
provide an ensemble average. Thus, it is critical to choose a task that may be made 
repeatably. After preliminary tests, it was found that a task that subjects could 
repeat accurately was to make rhythmic movements between two target points, in 
time to a regular auditory stimulus. This was the task performed in all experiments 
described. 

The two targets used were 1.0 rad (approximately 57 degrees) apart, well away 
from the joint limits (Figures 9-1 and 9-2). The upper arm was immobilized in the 
sagittal plane perpendicular to the torso. Target 1 was at the position with the 
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Target 1 
(0.0 rad) * 


0 



Target 2 
(1.0 rod) 


r (due to airjet) 


Figure 9-1: Target and variable assignments for movement task. 

forearm perpendicular to the upper arm, and target 2 was at 1.0 rad extension from 
target 1. The joint angle 9 is defined to be zero at target 1 and to increase with the 
elbow extension; thus, target 2 is at 9 = l.Orad. A single movement is defined to be 
one cycle from target 1 to target 2 and back. 

The auditory stimulus was provided by clicks from a piezoelectric buzzer, at in¬ 
tervals of either 750 ms or 1000 ms. Subjects were instructed to move continuously 
and repeatedly between targets, and be at a target at the time of a click. Continuous 
movement, not accuracy, was stressed. Subjects were instructed to relax and move 
naturally, in spite of the perturbations. 

In some trials the subjects were seated with their forearm moving in a vertical 
plane with respect to gravity. In others, subjects moved in the horizontal plane, with 
no gravity. 

In addition to the two speed conditions (750 ms and 1000 ms movements) and 
the two gravity conditions, in some trials the subjects were loaded with a damper 
(approximately 1 Nm/rad/s ) that acted about the elbow joint. The force-velocity 
specifications for the damper are shown in the Appendix. 

For each condition 15 blocks of 30 s each were collected (this gives 300 movements 
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Figure 9-2: Experimental setup. 

of 750 ms or 225 movements of 1000 ms per condition). Three minutes rest between 
blocks was provided to avoid fatigue. 

To control for joint angle dependent stiffness (e.g., due to muscle moment arm 
changes) subjects also maintained one of three poses while being perturbed (for 5 s 
each). These poses were at target 1, target 2 and a target midway between 1 and 2 
(i.e., at 0.5 rad). As in the movement experiments, subjects were asked to relax, but 
still achieve the task (i.e., maintain the target position). 

9.2 Stimulus Type: PRBS 

A pseudorandom binary force sequence, PRBS (see Figure 9-3), was chosen as the 
input stimulus because it could not be predicted by the subject, and it is optimal for 
system identification (in the sense that it can provide band limited white noise with 
a minimum length segment). See Chapter 8 for a review of the various alternative 
input stimuli. 

The PRBS was applied continuously by an airjet [150] [151] during all movements 
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Figure 9-3: Frequency response of a PRBS with a switching interval of 20 ms. 

tasks. The airjet applied forces of ±4 N to the wrist in the direction of elbow 
flexion/extension. 

In preliminary experiments the PRBS was composed of 10 ms elements (10 ms is 
the minimum stable switching interval for the airjet), giving a relatively white noise 
input signal up to 50 Hz. Subsequent analysis revealed that for this movement task 
the natural frequency of the arm was below 3 Hz (see results in next chapter). Thus, 
to put more low frequency power into the system at its natural frequency the airjet 
was switched at 20 ms intervals in a PRBS. This results in a relatively white noise 
input between 0 and 25 Hz. See Figure 9-3 for the frequency response of this PRBS. 
Frequency spectral shaping is possible by stochastic shuffling [154] or by analytic 
techniques [155], but was not deemed necessary without more a priori information on 
the elbow joint dynamics. 


9.3 Airjet Perturbation Device 

9.3.1 Specifications 

The primary design goal was that the PRBS perturbation device must not constrain 
the arm movement, or significantly change the arm dynamics. For this reason, a fight 
weight wrist-mounted pneumatic thruster device was used [150]. A further design 
constraint was that the device must have sufficient force and bandwidth to perturb the 
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arm measurably at frequencies up to 100 Hz. The 100 Hz frequency goal was chosen 
to be several times higher than the highest natural frequencies measured for joint 
dynamics (see measurements in [141] for example). This frequency allows accurate 
estimation of the inertia under all physiological conditions, provided the force is large 
enough. An approximately 4 N force turns out to be a good compromise for the the 
force produced by the airjet. As the inertia dominates at high frequencies it is possible 
to make a rough calculation of what displacement this force will produce at 100 Hz. 
With an inertia of I = 0.1 kgm 2 and a moment arm of 0.4 m, a 0.1 mm displacement 
results from a 4 AT 100 Hz force perturbation. This 0.1 mm displacement is within 
the resolution of the position sensor described below. Also, given the low frequency 
content of the PRBS perturbation stimulus it was found that a 4 AT force rarely 
perturbs the arm more than 5 — 10 degrees (see Chapter 10). Thus, the small 
perturbations assumption in the analysis is probably not violated. Finally, a 4 AT 
force applied at the wrist produces approximately the same torque as the maximum 
torque produced by the muscles in the task studied; see Section 11.1 for a discussion 
of this. Thus, the force is sufficiently large to provide relevant measurements to study 
normal speed movement formation. 

9.3.2 Design 

The details of the airjet design used in this thesis may be found in [150] [151]. A brief 
summary of the design is included here. 

The airjet is a device that is mounted on the subject’s wrist. It consists of a nozzle 
to speed up the air velocity to about 200 m/s , and a switching mechanism to divert 
the air flow down one of two thin walled brass tubes that bend to point in opposite 
directions along the line of arm movement. See Figure 9-4. With the flow down one 
tube a 4 AT reaction force is produced — due to the mass flow rate of air (30 cubic feet 
per minute). The force is switched 180 degrees to the other tube in 5 ms with a rise 
time of 1 ms. Figure 9-5 shows this switching response. The 5 ms delay is mostly due 
the solenoids used in the switching mechanism and could be lowered. The extremely 
rapid 1 ms rise time is realized because the switching is based on a fluidic property 
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inlet u, 


Figure 9-4: The wrist mounted airjet used in experiments. Figure reprinted from 
[ISO]. 



Time s 


Figure 9-5: Step response of airjet. Figure reprinted from [150]. 
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Figure 9-6: Force power spectrum of airjet following a PRBS at the maximum switch¬ 
ing rate. Note: this is not the force input used in the experiments. See Figure 9-3. 
Figure reprinted from [150]. 





sensor 


cuff 


Figure 9-7: Cuff mounting. The wrist is placed in the circular section, and the airjet 
attaches to the I-beam force sensor. Figure reprinted from [150]. 
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called the Coanda effect. Only the inertia of the air limits the theoretical switching 
speed. Figure 9-6 shows the device’s bandwidth. The mass of the airjet is 400 g. The 
air flow produces substantial noise, requiring the subject and experimenter to wear 
ear protection. 

9.3.3 Wrist-Cuff Attachment 

The airjet is placed above the hand such that it produces force in the direction 
of elbow flexion/extension and is aligned with the long axis of the arm to avoid 
pronation/supination torques. It is attached to an aluminum frame which is clamped 
on the wrist. See Figure 9-7. The subject wears a custom molded cuff for ease of 
clamping and comfort. Cuff movement was found to be negligible (0.1 - 0.2 mm) for 
the perturbations used [150], 


9.4 Sensing 

9.4.1 Force Sensing 

An I-beam with strain gauges connects the airjet to the cuff frame [150]. See Figure 
9-7. The I-beam is sufficiently compliant to obtain good force sensitivity, but stiff 
enough to have a resonance (160 Hz) well above the arm’s natural frequency. To 
avoid aliasing the force data were low pass filtered with a linear phase, —80 db stop- 
band analog filter with a cut-off frequency of 200 Hz (i.e., an 8 pole, 6 zero Frequency 
Devices DOW848 filter with a 4.2 ms pure delay). 

9.4.2 Position Sensing 

The position of an infrared light-emitting diode (IRED) on the cuff was measured with 
an Optotrak™ system, a 3D motion measurement system with a 0.05 mm resolution 
[150]. Data was sampled at the maximum rate of 200 Hz, and interpolated to the 
600 Hz force sampling rate with cubic splines. 

Let ( x , y) be any point in the plane of rotation of the elbow. The center of rotation 
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(x c , y c ) and radius r c to the IRED were estimated by finding the least squares best 
fit circle that goes through all the movement data. The angle with respect to vertical 
was then computed as: 


0 = tan ~ 1 ^—^ 
y~Vc 

The applied torque was computed as: 


(9.1) 


r = rF (9-2) 

where F is the measured force, and r = r c 4- 0.18m (0.18 m is the distance from the 
IRED to the force sensor). 

9.4.3 Pronation/Supination 

To control for possible unmeasured wrist pronation/supination some trials were per¬ 
formed with the airjet and cuff directly connected to a revolute joint with the axis 
of rotation aligned with the subject’s elbow. This eliminated possible wrist prona¬ 
tion/supination. In most subjects no difference was seen with and without the rev¬ 
olute joint support. Data with significant pronation/supination were thrown out. 
Ideally, additional IREDs should have been used to correct for non-planar motion, 
but the current Optotrak™’s limited bandwidth did not allow this. 

9.4.4 Data Acquisition and Filtering 

A real-time VME-based microprocessor system [156] was used to collect all data and 
also control the experimental stimuli. The sampled data were usually low pass filtered 
off line at 50 Hz. Additional filtering was applied as needed. All filtering was based 
on a recursive 12 th order Butterworth filter with a normalized cutoff frequency of 0.1 
[157, pg. 218]. Low pass, band pass and high pass filters of arbitrary frequencies were 
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obtained with the transformations in [157, pg. 238]. Phase distortion was eliminated 
by running the filter over the data twice: once in the normal manner, and then again 
with the data time reversed. 


9.5 Perturbation Calculation 

9.5.1 Alignment and Removing the Mean 

For each subject and movement condition, joint angle position profiles were obtained 
by (1) chopping the data into single ‘movements’ 0,(f) of one cycle each (where i — 1 
to n, and n is the number of movements made), (2) initially aligning the movements by 
peak velocity, (3) computing a mean movement profile, (4) re-aligning each movement 
so as to minimize the root mean square, RMS, difference between the mean profile 
and the movement 0j(i), and (5) re-computing a mean movement profile with the 
re-aligned data to give the final mean movement profile #o(t). 

Steps (4) and (5) were sometimes repeated to improve the average RMS differ¬ 
ence between the movements and the mean, but iteration was usually not neces¬ 
sary, and did not affect the final parameter estimates. Corresponding torque pro¬ 
files were obtained by chopping torque data up at the alignments computed for the 
angle data, computing a mean torque profile, and removing this mean from each 
torque profile. The final result gives n single movement position and torque tra¬ 
jectory profiles where i = 1 to to), a mean position profile and a mean 

torque profile (0o(t), To(t)), and to movement perturbation position and torque profiles 
(A0j(t), Arj(t), where i = 1 to to). Figure 10-2 shows typical movement data before 
and after removing the mean. 

9.5.2 Posture Perturbation Estimates 

In the case where the subject was asked to maintain a fixed posture, the dynamics 
were linearized about the mean arm position, 0 O . Thus, Ad(t) = 0{t) - 6 0 . Likewise, 
A r{t) = r(<) - r 0 . 
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9.6 Time-Varying System Identification 


The perturbation dynamics (7.9) were identified by finding the best parameter esti¬ 
mates J(t), B{t), and K(t) that satisfy the ensemble of measured perturbations at 
each fixed point in time. The method is summarized as follows: (1) (7.9) is dis¬ 
cretized^ j2) an optimization criteria is specified, (3) the discrete time parameter are 





9.6.1 Discretization 

The sampled data collected in the experiment should satisfy a general discrete time 
auto-regressive moving average (ARMA) model: 

0 = a o (t)A0(t) + ai (t)A8(t -h) + a 2 (t)A9(t - 2h)... 

+b 0 (t)AT(t) + & x (t)Ar(f - h) + b 2 (t)AT(t - 2 h) + ... (9.3) 


where h = 1.67ms is the sampling interval, the are the auto-regressive parameters, 
and the b{ are the moving average parameters. 

To relate the continuous time model (7.9) to (9.3) a simple method is to use the 
Euler approximation of a derivative: 

0(t + k)~ 0{t) 
h 

Substituting (9.4) into (7.9) gives: 




Ar(t) = a 0 (t)A9(t) + ai A0(t + h) + a 2 (t)A0(t + 2h) 


(9.5) 













































a sampled data zero-order-hold (ZOH) process, and integrate (7.9) over each sample 


interval [158], giving: 


0 = A 0(t — h) + ai(t)A0(t — h) + a 2 (t)A 6 (t — 2 h) 

-\-bi(t)Ar(t — h) 4 - b 2 (t)Ar(t — 2 h) (9-6) 

where the relation between ( I(t ), B(t), K(t)) and (ai(t), a 2 (<), &i(£), 62 ( 2 )) is non¬ 
linear (see page 52 in [158] for the inverse of the relation below): 

a i{t) + «2 (0 + 1 

W (t) + 

MW ( ,.t) 

Wo 

m 

Wq 

where z = + w 0 = w^/T+J, f = (1 fag - K{t)bx(t)lag - b/g ) 2 , g = 

sin(wh ), w = cos _ 1 ( 6 )/h, b — —ai(t)/ 2 a and a = ^Ja 2 (t). 

Exogenous variables [158] for modeling the noise process were not used. Instead, 
the following parameter identification method was found to be adequate. 

9.6.2 Parameter Identification with Input and Output Noise 

At each fixed point in time t the ARMA parameters («x(t), a 2 (t), &i(i), b 2 (t)) in (9.6) 
are constant for all n movements (where n = 300 for the 750 ms movements). These 
n linear equations in four unknowns may be solved in a number of ways. 

A standard method is to assume that all of the error enters additively into the 
output A 9(t), and find the parameters that minimize the mean squared error between 
the measured and predicted joint angles: 

E(A0<(f) - A 0^ l {t)f (9.8) 


K(t) = 

m = 
m = 
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where for the model (9.6): 


_A 0^\t) = a 1 (t)M i {t-h) + a 2 {t)M i {t-2h) 

+& 1 (t)Ar i (< — h) + b 2 (t)ATi(t — 2 h)) (9.9) 

and the subscripts i indicate the movement number. The problem with this linear 
least squares (LLS) method is that even if the measurements were corrupted only by 
ideal Gaussian white noise (GWN), the parameter estimates are biased because of 
(a) the auto-regressive terms (a x and a 2 ) and (b) input as well as output noise [158]. 

An improved method is developed by assuming that the input and output are cor¬ 
rupted by equal variance independent GWN [159]. Under these assumptions Koop- 
man [159] developed a simple unbiased estimator of the ARMA parameters (though 
not optimal in the Maximum-Likelihood sense). The method is analogous to the pro¬ 
cedure for fitting a multi-dimensional plane to data. Rewrite (9.6) for all movements 
as 


11 = 


A$ 0 (t) A 0 Q (t-h) A0 o (t-2h) Ar 0 (t — h) Ar 0 (t-2h) 
A $ 1 (t) A 0 1 (t-h) A0 1 (t-2h) Anit-h) Ar 1 (t-2h) 
• • a , t 
A0 n (t) A0 n (t-h) A0 n (t-2h) Ar n (t-h) Ar n (t — 2h) 


■ 

1 


«1 


d 2 


b x 


b 2 


or more compactly: 


(9.10) 


0 = Jp 


(9.11) 


Thus, p may be considered to be the normal to a plane in 5-dimensional space. Let 
the mean squared normal distance A from the data to the plane be: 
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A 


(9.12) 


= d 2 


p T J T Jp 

P r P 


The normal p that minimizes A is: 


p = the jth eigenvector of J T J 


(9.13) 


where the jth eigenvalue is the smallest eigenvalue (which also happens to be A, 
conveniently giving a measure of the model fit). The elements of p given by (9.13) 
provide the required ARMA parameter estimates. The eigenvector is readily obtained 
by singular value decomposition (SVD) of J without forming J r «7, thus avoiding 
numerical error accumulation associated with the J T J operation [58]. 

An estimate of the ratio of the input to output noise variance is necessary in order 
to scale the input measurements to achieve the assumed equal noise variance. For 
the data studied here, this ratio was typically set at 1.0, although changing it to any 
value between 0.01 and 100.0 did not appear to effect the results. 

The assumption that GWN is acting on the measurements is not as restrictive as 
it appears. As the measurements for each trial come from movements separated by 
more than one second it is reasonable to assume that noise across trials is independent. 
Only the data within three samples on a particular movement may be correlated (i.e., 
there may be correlated noise acting on A0<(f), A0j(f - h), A0,(t - 2h),Ar;(f - h), 
and A Ti(t — 2h)). 


9.6.3 Time-Varying Impulse and Frequency Responses 

The time-varying impulse response may be modeled using (9.3) with a 0 (t ) = —1, 
ai(t) = 0 for l > 0, and fy(i) not necessarily zero for lags l up to several times the 
system response time. Parameter estimation may be carried out by an appropriately 
modified (9.10), or by the linear least squares (LLS) method. Here LLS is chosen as 
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there are no auto-regressive terms. In matrix form (9.3) is: 


A0 O (t) 


A#i(t) 

= 

• 



At 0 (<) 

Arx(t) 


Ar 0 (t - h ) 
Ar 1 (< — h ) 


A r 0 (t — 2 h) 
A T X (t — 2 h) 


or more compactly: 


bo(t) 


(9.14) 


y = 


(9.15) 


and the LLS solution is 


Pi = (JfJi) 1 Jfy (9-16) 

Notice that JfJi and Jfy_ are respectively the ensemble auto- and cross-correlations. 
The time-varying frequency response can be estimated by Fourier transforming the 
impulse response at each fixed point in time. 

The problem associated with this analysis is that at each time instant there are at 
least 180 parameters (i.e., bi(t) with Z = 0 to 179 for a 300 ms long impulse response 
at the 600 Hz sampling rate), but only 200 to 300 equations (movements) to use in 
the estimation. Thus, as a result of noise, «7j typically is not of full rank and the 
solution (9.16) is undefined. 

To generate more equations it is reasonable to assume that the dynamics are 
time-invariant within fixed time windows w (w = 60ms). Thus, within a window n 
equations of the form (9.14) hold at each of the wf 3 (=37) sample times. Thus, there 
are wf t n (= 37n) equations with which to estimate the fixed ARMA parameters in 
a given window. 

Even with the windowing technique described in the previous paragraph, it was 
found that for the data obtained in these experiments J\ was still not of full rank. 
To uniquely define a solution the squared Euclidean norm of ^ may be minimized 
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simultaneously to minimizing the mean squared modeling error: 


(y - JiPifiv ~ JipJ + Apfp, (9.17) 

This minimization may be implemented by ridge regression [160], or equivalently 
with the singular value decomposition SVD approach described in [58]. (For the SVD 
method, any singular value less than 10 -4 of the maximum singular value was zeroed.) 
As this solution is only of mathematical convenience and the real problem is one of 
having too many parameters to be estimated, the impulse response estimates are of 
questionable value and were only estimated for a few data sets (see Figure 10-7 in 
the results section). Perhaps a better solution would be to minimize the norm of the 
discrete first or second derivative of the impulse response p^ (rather than the norm of 
£j), thus enforcing a smoothness constraint on the impulse response. 

9.6.4 Posture Impedance Estimates and Ergodicity 

In the case of the perturbation data collected during posture, the above system iden¬ 
tification technique must be modified as there is only one input/output pair. To do 
this it is assumed that the system parameters I(t), B(t), and K(t) are constants J, 
B , and K respectively. An equation of the form (9.3) with constant parameters thus 
describes the data. Therefore, combining the equations across time, rather than down 
the movement ensemble, is valid. With A0j(t) = A 6{t — i ) and A Ti{t) = A r(t — i) 
substituted into (9.10) the constant vector p may be estimated and /, B and K 
recovered using (9.7). 

Inspecting the elements of J T J reveals time-averaged cross-correlations and au¬ 
tocorrelations of the input and output variables. Effectively, the substitution of 
A $i(t) = A 6(t — i ) and Ar^(t) = Ar(t — i ) turns the ensemble averaged correlations 
into time-averaged correlations. Thus, in a statistical framework the justification 
for this substitution is ergodicity (i.e., time averages equal ensemble averages). No 
assumption about the stationarity of the input signal A r{t) is used, although its sta- 
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tionarity would make the autocorrelation matrix (a submatrix of J T J ) Toeplitz in 
structure, which could be exploited to simplify the computation. 

Analogously to the previous section, the impulse response during posture may be 
modeled by assuming a time-invariant system with a structure of (9.3) with a 0 — — 1, 
ai = 0 for l > 0, and bi not necessarily zero for lags l up to several times the system 
response time (typically, at least 300 ms, or l = 0 to 179 for the 600 Hz sampling 
rate). Parameter estimation may be carried out using an appropriately modified 
(9.10). Alternatively, as no autoregressive terms exist to introduce parameter bias, 
linear least squares LLS may also be used (assuming no input noise). The LLS 
solution is equivalent to the standard correlation solution. The frequency response 
may be estimated by Fourier transforming the estimated impulse response. 

9.6.5 Simulations 

A critical part of any system identification technique is verification with simulated 
data. To simulate the experimental setup the forearm was modeled as a cylinder (mass 
and center of mass m and r) with inertia I = 0.07 bkgm 2 , the muscles were modeled as 
linear with constant stiffness K = 20Nm/rad and damping B — 1.0 Nm/rad/s, and 
the arm was considered to move in the vertical plane against gravity. The airjet was 
assumed to have a mass of m a = 0.2kg acting at r a = .4m radius, giving an inertia 
of I a = ,002kgm 2 . Thus, without the airjet running the unperturbed dynamics are: 


0 = I0(t) + IJ(t) + B0(t) + K(0(t)-0 o (t)) (9.18) 

—mgr sin 0{t) — m a gr a sin 0(t) (9.19) 

The motor input 0o(t) is calculated so that the unperturbed dynamics (9.19) pro¬ 
duce an arm trajectory 0(t) that is a 0.5 Hz sinusoid with a peak-to-peak amplitude 
of 1 rad (similar to the 1000 ms movement condition, where each cycle is one move¬ 
ment), and 100 cycles long. 
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With the airjet running, the perturbed dynamics are: 


r(£) = I0(t) + B0(t) + K(0(t) — 0 o (t))—mgr sin 0(t) (9.20) 

where r(t) = rpRBs(t) — l a 0{t) — m a gr a sin 0(t) is the force measured by the force 
sensor on the airjet mounting, and TpRBs{t) is a PRBS of ±1.6 Nm. Using 0 O (t) and 
r PRBs(t ) as input to a fourth order Runge-Kutta numerical integration [58] of (9.20), 
the perturbed joint trajectory 0{t) and perturbations measured by the force sensor 
r(t ) are computed. 

The above procedure was used to generate simulated data. The data was then 
processed as in the experiments. That is, the simulated data 0{t) and r(£) were 
chopped into single movements (one cycle each), and aligned with the technique 
described in Section 9.5.1. Mean movement profiles r 0 (£) and 0o(t) and perturbations 
about the mean were computed. See Figures 9-8 (a), (b) and (c). 

The perturbation parameters were then estimated using (9.13) (after introducing 
5% GWN noise into the input and output data). See Figure 9-9. Notice that only the 
inertia of the arm I was estimated, not the total inertia of the airjet and arm. Also 
notice that the stiffness estimates are a little lower than K = 20 Nm/rad, particularly 
near target 1 (0 = 0). This is the effect of the negative gravity contribution to stiffness, 
as discussed after equation (7.9). 

9.6.6 Apparatus Identification 

The mass of the wrist cuff and the frame supporting the airjet contributes to the 
measured inertia. To estimate the inertial contribution of the whole apparatus it 
was supported on a light low friction revolute joint, such that the airjet was at the 
same distance and orientation relative to the joint that it typically is during the 
experiments. A PRBS input was applied by the airjet and measured as in Section 
9.4. The acceleration was computed by 50 Hz low pass filtering the recorded joint 
angles and using the approximation (9.4) twice. The static relation r = 1^/0 was 
then fit to the data to provide an apparatus inertia estimate of /«*// = 0 . 035 &< 7 to 2 . 
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9.6.7 Artificial Arm Identification 


Using the same support joint for the airjet apparatus described in the previous section 
a simulated arm maybe made by coupling the joint to a damper and/or spring. Using 
only a damper (see Appendix) a simulated experiment was performed where the airjet 
was switched with a PRBS with 20 ms elements. The method of section 9.7.2 was 
used with (9.13) to estimate the perturbation parameters, giving: I = 0.04fc^m 2 , 
B = 1.20Nm/rad/s and K = 0.91 Nm/rad. These estimates agree well with the 
inertia estimate of the previous section and the specification sheet provided by the 
manufacturer of the damper (see Appendix). In an independent test a 4 N step 
response also gave a similar damping estimate (1.3 Nm/rad/s). The small stiffness 
could be due to the restoring forces provided by the air supply lines. 
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Chapter 10 


Results of Airjet Perturbation 
Experiments 

10.1 Posture 

The angle dependence of perturbation parameters estimated during posture (see Sec¬ 
tion 9.6.4) is shown in Figure 10-1. The plotted values show the mean and standard 
deviations of six estimates made from separate 5 s trails where subjects were in¬ 
structed to point at a fixed target location. For the horizontal movements there is 
little dependence of these perturbation parameters on the joint angle. For the vertical 
movements the stiffness increases slightly with joint angle. 

The torque variance accounted for (VAF) by the second order perturbation model 
was computed for each trial, and is consistently above 80% (where the VAR/100% is 
defined as one minus the ratio of the variance of the model’s error in estimating torque 
to the variance of the torque data). Frequency responses and coherence functions were 
calculated for all subjects. The coherence was always above 95% up to 50 Hz. 

10.2 Horizontal 750ms Movements 

In this and following sections we now consider how the mechanical properties are 
modulated under various movement conditions. First we focus on the horizontal 
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Figure 10-1: Joint angle versus mean (a) inertia, (b) viscosity and (c) stiffness for 
Subject SI asked to maintain fixed postures in the vertical plane (marked with circles) 
and horizontal plane (marked with triangles), (d) Joint angle versus mean stiffness 
for Subject S2. Error bars show standard deviations about the means. 
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plane 750 ms movements. 


10.2.1 Typical Movement Data 

Typical angle and torque measurements are shown in Figure 10-2 for two subjects. 
Figure 10-3 shows the corresponding perturbations computed by removing the en¬ 
semble means (re Section 9.5.1). 

10.2.2 Main Results 

The time-varying parameter estimates made using the perturbations from 300 move¬ 
ments and (9.13) are shown in Figure 10-4 (See Figure 10-6 also). The mean parameter 
values and standard deviations computed over 60 ms windows are shown. The data 
are repeated for an extra cycle to emphasize that the estimates are periodic. The 
corresponding damping ratios (B/(AKI) 0 - 5 ) are shown in Figure 10-5. The features 
in Figures 10-4 and 10-5 that are typical across all subjects and conditions are: 

1. The stiffness is modulated with the movement. It reaches a peak value just before 
the arm comes to rest at a target and drops to a minimum value before the peak 
velocity of movement is reached. 

2. The stiffnesses are low, lower than the smallest values measured during posture. 

3. The stiffness values at the two targets are sometimes unequal (usually higher at 
target 2 when gravity is present). 

4. There is a tendency for the damping to drop slightly while the arm is at a target. 
The damping estimates have considerable variation across subjects though. 

5. The damping ratio is not constant, ranging from 0.2 to 0.6 during the movement. 

6 . The inertia estimates are constant and agree with previous estimates of forearm 
inertia. The values shown include the apparatus inertia of O.OSSkgm 2 . 
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Figure 10-2: Typical horizontal plane elbow angle and torque profiles for subject SI 
and subject S2. Each plot represents a full cycle taken from a rhythmic oscillation of 
the forearm in time to a periodic auditory stimulus (period 750 ms and duration 30 
s). Targets to be hit were at 0.0 rad and 1.0 rad. 


L48 



Figure 10-3: Typical horizontal plane elbow joint perturbations about the ensemble 
means for subject Si and subject S2. Perturbations were computed by subtracting 
the mean of 300 movements from the trajectories in Figure 10-2. 
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Figure 10-4: The estimated time-varying perturbation model for subject SI moving 
in the horizontal plane in time to a 750 ms auditory stimulus. The mean parameter 
values and standard deviations computed over 60 ms windows are shown. The data 
are repeated for an extra cycle to emphasize that the estimates are periodic. 
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Figure 10-5: Figure 10-4 continued. Damping ratio profile. 

The stiffness variations mentioned above cannot be simply due to passive connec¬ 
tive tissue or skeletal moment arm changes, as the stiffness values measured during 
posture do not vary with joint angle (over the range studied). 


10.2.3 VAF 

The torque variance accounted for (VAF) by the model parameters in Figure 10-4 is 
76%. Across all conditions and subjects the VAF by the time-varying second order 
system parameters exceeds 70%. Most of this error is due to voluntary variation in the 
trajectory from movement to movement, which effectively introduces low frequency 
noise into the data. 

10.2.4 Frequency Responses and Natural Frequencies 

Three samples of the time-varying frequency response (at (a) 9 = 0.0 rad, (b) 9 — 
0.5 rad on the extension movement, and (c) 9 = 1.0 rad) are shown in Figure 10-7. 
These were estimated with the SVD method of Section 9.6.3 with a 60 ms window, 
and a Fourier transformation. The natural frequencies, as defined by the point at 
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which the phase angle crosses -90 degrees, are respectively (a) 1.77 Hz , (b) 0.63 Hz, 
and (c) 2.51 Hz. The corresponding natural frequencies computed from the second 
order system parameters are respectively (a) 1.50 Hz, (b) 1.06 Hz, and (c) 2.05 Hz 
(where the natural frequency f n is defined as 2i r/ n = (if//) 0,6 ). This second method 
of natural frequency estimation is superior in view of the problems mentioned earlier 
in connection with estimating the time-varying impulse response. By all measures, 
the natural frequency is always below 3 Hz. 

10.2.5 Dominance of Inertia Above 5 Hz 

To demonstrate the dominance of inertia above 5 Hz the data were filtered with a 
high pass filter with a 5 Hz cutoff frequency (see Section 9.4.4), and a static inertia 
model Ar = IAS was estimated. Figure 10-8 shows the measured torque (thin line) 
and the estimated torque I Ad (dotted line) for a segment of the data. Both angle and 
torque data were also low pass filtered at 30 Hz in order to differentiate the angle data 
to get acceleration — thus the smoothness of the curves. The variance accounted for 
(VAF) by this constant inertia model was 95%, verifying that the system’s natural 
frequency is below 5 Hz, and indicating that the inertia dominates the dynamics 
above 5 Hz. 


10.2.6 Effect of Airjet Perturbation Frequency- 

Figure 10-9 shows results of an experiment in which subject S2 moved in the horizontal 
plane and was perturbed with the airjet switched at twice the usual frequency (10 
ms elements in the PRBS). With this stimulus rate the input was white up to 50 
Hz, but there was only half the power that there usually is between 0 and 25 Hz. 
Comparing Figures 10-6 and 10-9 shows that this increase in bandwidth and decrease 
in low frequency amplitude does not qualitatively affect the results. Quantitatively, 
the VAF in Figure 10-6 is 72% and the VAF in Figure 10-9 is 76%. 


153 



m 



Mlll l UlM I.IMM 

#*_»/** (Ml) 


Figure 10-7: Time-varying frequency responses at (a) 9 = O.Orad, (target 1), (b) 
9 = 0.5 rad (on the extension movement), and (c) 9 = 1.0 rad (target 2). The natural 
frequencies, given by the point at which phase angle crosses -90 degrees, are respec¬ 
tively (a) 1.77 Hz , (b) 0.63 Hz, and (c) 2.51 Hz. Data are for subject S2 moving at 
the 750 ms rate in the vertical plane. 
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Figure 10-8: The data from Subject Si’s horizontal movements (750 ms stimulus) 
were filtered at 5 Hz with a high pass filter, and a static inertia model At = IA6 was 
estimated. Plotted above are the torques (thin line) and the estimated torques I AO 
(dotted line) for a segment of the data. The variance accounted for (VAF) by this 
constant inertia model is 95%, indicating that the inertia dominates the dynamics 
above 5 Hz. 

10.3 Effect of Speed Change 

The stiffness estimates increase slightly with a speed increase. Figure 10-10 shows 
typical results when moving between targets in the horizontal plane at the 0.5 Hz rate 
(1000 ms stimulus). Comparing this to the equivalent but faster 0.67 Hz movements 
in Figure 10-4 it is seen that the stiffness profiles are similar. The stiffness does not 
scale up by a factor equal to the square of the speed-up factor, which is x 1.78. Figure 
10-11 (also see Figures 10-15 and 10-17) compares the stiffness profiles for the two 
speeds for subjects SI and S2. 


10.4 Effect of a Viscous Damping Load 

Figure 10-12 shows the results of attaching a damper to the arm during movement. 
Comparing these results to Figure 10-4 reveals that the damper more than doubles 
the damping that is normally measured in the unloaded joint (note that the damper 
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Figure 10-10: Same as Figure 10-4, except with a 1000 ms period auditory stimulus. 
Data are for Subject Si. 
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has a non-linear force-velocity relation shown in the Appendix — thus the higher 
damping at the targets). In spite of this increase in damping, the peaks in the 
stiffness profile are only increased by 20 percent, and the minimum stiffness during 
movement is unchanged. Effectively, the system’s damping ratio is doubled, and there 
is no attempt to compensate. 



Figure 10-12: Same as Figure 10-4, except with a viscous damper loading the arm. 
Notice that the scale on the damping is different. The damper load is non-linear 
(the damping drops slightly with increased velocity. See Appendix). Note that the 
stiffness profile is very similar to Figure 10-4. The inertia is a little larger than in 
Figure 10-4 because of the damper attachment. Data are for Subject SI. 
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10.5 Effect of Gravity 

10.5.1 Stiffness Increases Due to Gravity 

To observe the affect of gravity, analogous experiments to the horizontal plane ex¬ 
periments just described were performed in the vertical plane. Figures 10-13 and 
10-14 show the results from subjects Si and S2 moving in the vertical plane in time 
to the 750 ms auditory stimulus. Comparing these to the horizontal movements in 
Figures 10-4 and 10-6 reveals that the stiffness is higher at target 2 in the vertical 
plane movements. Note that the inertia estimates in the vertical movements are still 
constant, but are slightly higher in both subjects. This difference can be attributed to 
the cuff being mounted slightly further down the arm. This shift in mounting should 
not affect the stiffness or damping estimates. 

To see the stiffness profile differences more clearly, the stiffness profiles from two 
subjects and speed conditions were plotted together in Figure 10-15. Comparing 
these results to the corresponding horizontal plane plots in Figure 10-11 reveals that 
there is a significant asymmetry in the vertical plane stiffness profiles. Unlike the 
horizontal plane movements the stiffness at target 2 (1.0 rad ) is two to three times 
higher than at target 1 (0.0 rad). Notice that at target 2 the subject must resist the 
larger gravity force. This will be elaborated further in the next chapter. 

As mentioned in Section 7.3 there is a small contribution to the measured stiffness 
that is due to gravity (—mgc cos 9). At vertical {6 = O.Orad) this contribution is 
highest. To see the size of this contribution it was computed for each angle and 
subtracted from the measured joint stiffnesses. The required mass and center of mass 
of the arm were estimated from the measured inertia using an assumption that the 
arm could be modeled by a cylinder. Figure 10-16 shows the stiffness profiles (a) 
before and (b) after making this gravity correction. Note that the difference does not 
account for the asymmetry in the stiffness peaks. As the difference between these two 
profiles is small, this gravity contribution will not be considered further. 
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Figure 10-14: Same as Figure 10-6, except that the movement occurred in the vertical 
plane. Data are for Subject S2. 
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Figure 10-16: The negative gravity contribution to stiffness is small. Plot (a) shows 
the measured stiffness profile for the subject moving in the vertical plane and (b) 
shows the same profile with the gravity contribution to stiffness (—mgc cos 8) removed. 
Data are for Subject S2. 
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10.5.2 Effect of Speed Change Under Gravity 

As with the horizontal movements, there is only a marginal increase in stiffness with 
an increased speed of movement in the vertical plane. Figure 10-17 shows the stiffness 
profiles for four subjects moving in the vertical plane at the two speeds. 



Figure 10-17: Stiffness profiles measured in the vertical plane at the two speeds. 
Subjects as indicated. 
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Chapter 11 


Discussion 


11.1 Relevance of Perturbation Estimates 

An important issue that must be cleared up before proceeding is whether or not 
the mechanical impedance estimates are relevant to the formation of normal speed 
movement. 


11.1.1 Short Range Stiffness 

It might be argued that the low amplitude, high frequency input of the airjet only 
excites higher order dynamics that are not normally relevant to unconstrained move¬ 
ment (e.g., short range stiffness [105]). This argument is wrong for several reasons. 

First, care was taken to use a sufficiently large perturbation force. In fact, the 
1.6 Nm torque produced by the 4 N airjet thrust is comparable to the maximum net 
torque that the muscles produce in the movements studied, see Figure 11-1. Thus, for 
the segments of the PRBS in which the airjet is on in one direction for an appreciable 
portion of the movement time (say >100 ms) the arm effectively receives a step torque 
change comparable to the maximum net muscle torque. 

Second, after the initial experiments that indicated that the arm’s natural fre¬ 
quency was below 3 Hz , the airjet’s minimum switching time was raised to 20 ms, 
providing more power at the low frequencies. Figure 9-3 shows that the majority of 
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the input power is evenly distributed between 0 and 25 Hz. As a qualitative indica¬ 
tion of the low frequency power, it is mentioned that the airjet regularly perturbed 
the arm 5 to 10 degrees, and on occasion, for the arms of smaller subjects, displaced 
the arm the entire 60 degrees between targets (in these few cases the subject would 
abandon that movement). 

Finally, the stiffnesses and corresponding natural frequencies measured in this 
thesis are much too low to be of only short range stiffness origin. As mentioned, 
the stiffness estimates give values for the system’s natural frequency below 3 Hz. A 
separate corroboration of this low natural frequency comes from noticing that 95% 
of the variance may be accounted for by a pure mass model, if the input and output 
are high pass filtered above 5 Hz (Figure 10-8). If the airjet was probing only the 
short range stiffness, much higher natural frequencies estimates (stiffness estimates) 
would be expected. For example, the short range forearm stiffness estimates made 
by Lanman [130] (with high frequency 20-30 Hz sinusoidal inputs) are 5 to 10 times 
the values estimated in this thesis. 

11.1.2 Spindle Reflex Sensitivity with Vibrations 

The presence of the low amplitude high frequency content in the input PRBS may al¬ 
ter the muscle spindle sensitivity. The spindle response may be saturated (or clamped) 
by the vibrations, making them insensitive to larger, lower frequency stretch [161]. 
Perhaps there is a connection between this and the observation made in Chapter 8 
that the damping ratios measured with PRBS inputs are lower (0.3 to 0.6) than those 
measured with single pulse inputs (0.8 to 1.0). There is no simple method to control 
for this problem, as the apparatus used can only produce binary force sequences. As 
much of the following analysis relates to relative stiffness changes with movement, 
these possible spindle sensitivity changes do not affect the conclusions. 
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11.2 Testing of Control Theories 


The predictions (Pi) to (P5) of Chapter 7 are now tested with the perturbation 
parameter measurements from Chapter 10. 

11.2.1 Compliant Feedforward Control 

First consider prediction (PI). It is clear from the horizontal movements in Figure 
10-4 (and others like it) that, as predicted, the joint stiffness is low and only rises 
as the targets are approached. As discussed earlier, the finding that the stiffness is 
low suggests that feedforward control dominates. Further, joint stiffness rises exactly 
when it is required to reject disturbances in attaining the target; that is, it reaches 
a peak just before the target is attained. Once at the target no further accuracy is 
needed and the stiffness can drop immediately. 

11.2.2 Stiffness and Force Modulation 

As the arguments behind the compliant feedforward control hypothesis are purely 
functional, it is interesting to inquire into the mechanism behind the stiffness mod¬ 
ulation. This will be attempted in the next section, but to start toward this goal 
we make the hypothesis that the stiffness changes are approximately proportional to 
muscle forces changes. 

First, consider the horizontal movements. Assuming little co-contraction, the 
torque from the agonist muscles can be calculated to be the torque needed to accel¬ 
erate the arm. As the acceleration is a maximum at the targets, the muscle torque is 
also at a maximum. Figure ll-l(a) shows the muscle torque required to accelerate the 
arm inertia (computed from the estimated inertia and average acceleration). Notice 
that indeed the stiffness goes up where the muscle torque is highest. 

Further, comparing the two speed movements it is found that the stiffness is only 
higher in the faster movements near the targets. This is also true for the comparison 
between the computed net muscle torques at the two speeds; the difference is largest 
near the targets. In the middle of the movement the acceleration is zero, irrespective 
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of the speed. See Figure ll-l(a-b). The stiffness increase is not sufficient to conclude 
a strict proportionality with force, nevertheless the stiffness does go up with muscle 
force increase. 

Likewise, the muscle torque in the vertical movements can be calculated. In this 
case gravity also contributes to the torque. Figure ll-l(c) shows the computed net 
muscle torque (where the computation uses the estimated inertia, a cylinder model 
of the arm to estimate the mass an center of mass from inertia, and the average 
trajectory). Notice that the stiffness is again highest where the muscle torques are 
highest. The muscle torques at the two targets are asymmetric, just as the stiffnesses 
are. 

Unfortunately, this simple explanation for muscle stiffness changes is not complete. 
Consider the damper experiment. This experiment was specifically designed to test 
this force-stiffness explanation. The damper produces its largest torque in the middle 
of the movement (at the peak velocity). Figure ll-l(d) shows the computed net 
muscle torques to oppose inertial acceleration and the damper load. But Figure 10- 
12 shows that the stiffness does not go up during movement, despite the large increase 
in torque during the movement. The stiffness does go up slightly, but only near the 
targets. Thus, the mechanism for stiffness modulation must be more complex than 
this simple monotonic force-stiffness relation. This may have been expected, as it 
was mentioned earlier that muscle force is velocity dependent. Also, reflex induced 
stiffness changes and co-contraction were ignored in the above arguments. 

11.2.3 Speed Scaling and Feedforward Control 

To further test the hypothesis that feedforward control dominates, consider the dif¬ 
ferent speed movements again. If there was no feedforward compensation (other than 
perhaps static feedforward gravity compensation) then the effective feedback gains 
(as measured by the stiffness) would have to increase with the square of the speed 
increase (i.e., xl.78) to produce the same time scaled movement trajectory (i.e., pre¬ 
diction (P2)). In actuality, time scaled movement trajectories at the two speeds are 
the same (e.g., Figure ll-2(b-c)), and yet for the faster movement the peak stiffness 
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Figure 11-1: Torques that the muscles/tendons produce in subject SI during (a) an 
average horizontal 750 ms movements, (b) an average horizontal 1000 ms movements, 
(c) an average vertical 1000 ms movement, and (d) an average horizontal 750 ms 
movement with a damping load. 
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Figure 11-2: Mean trajectories (position and velocities) for subject SI during the 
(a) horizontal 750 ms movements, (b) vertical 750 ms movements, (c) vertical 1000 
ms movements, and (d) horizontal 750 ms movements with a damping load. All 
amplitudes are normalized. For reference, the 1000 ms and 750 ms peak velocities 
are respectively 1.6 rad/s and 2.0 rad/s , and the peak to peak displacement are 
approximately 1 rad (see previous figures). 
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is on average only a factor of 1.18 larger, and the mean stiffness across the movement 
is on average only a factor of 1.05 larger. These stiffness changes are significantly 
different from the predicted factor of 1.78. Thus, prediction (P2) is incorrect and a 
feedforward speed compensation is necessary. 

The dynamic feedforward compensation is not complete though. As predicted in 
(P4) the trajectories loaded with a 1 Nm/rad/s damper are different from the tra¬ 
jectories in all other conditions. Figure 11-2 shows this difference and also shows the 
similarity between the trajectories of all other conditions (gravity and speed). Partic¬ 
ularly notice how different the velocity profile is for the damping loaded movement. 
It is perhaps not surprising that the arm does not maintain the same trajectory, as 
the task only requires reaching the targets. 

Finally, as stiffness and viscosity do not scale up appropriately with speed, the 
possibility that the feedforward compensation is simplified by dynamic scaling is 
disproved. That is, prediction (P5) is incorrect. The stiffness and damping are not 
scaled so as to allow simple scaling of the feedforward torque command to produce 
faster movements. 


11.2.4 Static Feedforward Gravity Compensation 

Though the feedforward compensation is not merely static, it is still possible that 
gravity compensation occurs, giving the prediction (P3). Comparing the vertical 
movements to the horizontal movements (e.g., Figures 11-2, 10-15 and 10-11) it is 
found that there is little difference in either of the trajectories and the only difference 
in the dynamics is a less than 5 Nm/rad mean stiffness increase. Thus, unless this 
stiffness difference is significant, it is probable that gravity compensation is being 
used. 

To further investigate the affect of these stiffness differences between the vertical 
and horizontal movements consider the arm to be driven by a linearized model of the 
form (7.15). In the absence of gravity compensation the computed desired trajectory 
0 O is shown in Figure ll-3(a). This was computed from (7.15) using the average actual 
trajectory, the estimated parameters, and a cylinder model of the arm (to estimate 
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the mass and center of mass from the inertia). This is the trajectory that the the 
equilibrium trajectory hypothesis [136] [28] predicts. With gravity compensation the 
effect of the last term in (7.15) is zeroed, and the computed desired trajectory is shown 
in Figure ll-3(b). Notice that without gravity compensation the desired trajectory 
must take on a complex shape. Of course, these plots must be treated with care 
because they presume linearity of the arm dynamics, or at least that the identified 
parameters approximate the non-linear dynamics adequately (see Section 7.4). In 
any event, it appears that feedforward gravity compensation is required. That is, the 
stiffness increase in the vertical movements is not enough to compensate for gravity 
without explicitly providing a feedforward term to cancel out the static effect of 
gravity. 

11.3 Mechanisms for Stiffness Modulation 

As mentioned before, the perturbation stiffness measurements are a function of passive 
joint properties, muscle, and reflexes. Though it is difficult to distinguish these factors 
there are certain muscle properties that may well aid in explaining the changes in 
stiffness with velocity. 

11.3.1 Intrinsic Velocity Dependent Compliance 

Consider the repeated horizontal movements studied here. In this task the arm ex¬ 
ecutes roughly a sinusoidal 0.5-0.7 Hz trajectory by exerting reciprocal biceps and 
triceps activity. A given muscle, when activated, can be considered to be driven by 
a sinusoidal position drive (i.e., the drive resulting from arm inertia and the other 
muscles). Thus, the isolated whole muscle oscillation studies of Rack and Westbury 
[105] are relevant to the present study. They oscillated tetanized cat soleus muscle 
at 0.9 Hz, and plotted the instantaneous length versus tension curves for complete 
cycles of the input sinusoid (the cycles were non-conservative: encircling a non-zero 
area in the length-tension plane, requiring work). Starting from zero velocity, the 
initial stretching of the muscle produces a sharp increase in tension until the stretch 
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Figure 11-3: (a) Actual trajectory (solid line) and desired trajectory computed under 
the equilibrium trajectory hypothesis (dotted line), (b) Same as in (a) but gravity 
compensation is assumed in addition to the equilibrium trajectory hypothesis. Data 
shown are for subject S2. Desired trajectories were computed with equation (7.15) 
as described in the text. 
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exceeds 1 mm, and then the muscle continues to stretch with little increase in tension. 
On the reverse part of the cycle exactly the same thing happens: starting from zero 
velocity, the initial shortening produces a symmetrical sharp drop in tension until the 
shortening exceeds 1 mm, and then the muscle shortens with little decrease in ten¬ 
sion. In effect, the muscle has a short range high stiffness acting only when starting 
from zero velocity, and only lasting for 1 mm stretches (5 percent of muscle length, 
or 4° of ankle rotation [122]). In the rest of the cycle (that is, during movement) the 
muscle is very compliant, resisting length changes with little force. 

This phenomena was also observed in their earlier studies [121] [122] [123]. Figure 
7-8 (the solid lines) shows their results for ramp stretches starting from various lengths 
and tensions. The ramp inputs cause an initial sharp rise in tension, which lasts only 
for about 1 mm (4° of ankle rotation). The tension then only rises slowly, or in some 
cases drops off, as the lengthening continues. 

It can thus be speculated that the airjet perturbations are met with high short 
range stiffness only when the arm is near rest, that is, within a few degrees (« ±4°) 
of the target. This short range stiffness is only probed by segments of the PRBS that 
switch fast enough to maintain the arm displacements within a few degrees. The low 
frequency segments of the PRBS perturb the arm outside of this short range stiffness 
zone, and would thus tend to lower the stiffness estimates. 

Further, it may be speculated that during movement the airjet PRBS induced 
displacements should meet with much more muscle compliance than when at rest. 
That is, the tension from shortening active agonist muscles should drop relatively 
little in response to the displacements in the direction of movement (i.e., additional 
small displacements in the direction of the already shortening agonist should affect 
the muscle force similarly to the plots in [121][122] [123]). Likewise, the tension in the 
lengthening antagonist muscles should rise relatively little in response to perturbation 
displacements in the direction of movement. Figure 7-8 indicates that it might even 
be possible for the tension to drop (for a lengthening muscle just beginning to move) 
in response to stretch. The muscle responses to PRBS induced displacements in the 
opposite direction to movement are probably similar if the displacements are small 
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enough not to reverse the direction of movement. 

Lanman [130] (see review in Chapter 8) probed the human arm with high fre¬ 
quency sinusoidal perturbations during movement, and also found that the muscle 
stiffness drops during movement (the frequencies he used are too high to allow reflex 
activity). He attempted to link these results to a molecular level description of muscle, 
arguing that perhaps for displacements above a certain threshold the actin-myosin 
bond breaking rate is increased in proportion to the velocity. Thus, as the velocity 
increases the number of bonds should decrease, and the stiffness provided by the sum 
of all bonds should decrease. 

Lanman also demonstrated that if the arm is loaded with a viscous damper during 
movement, then the stiffness no longer dropped. He likewise explained this result with 
a crossbridge model. He argued that the increased load requires more crossbridges 
to be formed, offsetting the effect of the crossbridges being broken by the muscle 
length changing. Lanman’s result (and explanation) is hard to reconcile with the 
experiments of [105] [122] described above, but is probably related to the fact that 
he only probed the arm with single high frequency (15-30 Hz) sinusoids. Also, as 
mentioned in Chapter 7 the distributed sarcomere force generation mechanisms are 
not generalized easily to whole muscle. 

Generalizing Lanman’s explanation to explain the data collected with PRBS per¬ 
turbations with significant low frequency content must proceed with caution. As 
mentioned at the beginning of this chapter, below 5 Hz the airjet produces pertur¬ 
bation length changes beyond the short range stiffness of muscle (e.g., short range 
stiffness is attributed to the elastic limit of crossbridges in [130]). Although the mea¬ 
sured stiffnesses do tend to drop with velocity, they do not increase much with a 
damping load. Thus, there must be other mechanisms than short range stiffness that 
maintain the stiffness profiles across these experimental conditions. 

In summary, known muscle properties could qualitatively account for some of the 
features of the data collected in this thesis. Muscle can be viewed to have long time 
constant (> 50 ms) elastic properties (re length-tension curves) if allowed to come 
to rest (see review of muscle properties in Chapter 7). Further, at rest muscle can 
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provide pure high stiffness if not stretched beyond the limit of its short range stiffness 
(e.g., elastic limit of the crossbridges). Both of these stiffnesses contribute to the 
measured stiffness near the targets. Also, as more force is required in the vicinity of 
the target (e.g., due to acceleration gravity changes) more parallel force generating 
elements have to be recruited (either by increased firing frequency or increased muscle 
fiber recruitment), and thus the stiffness should go up. On the other hand, during 
movement muscle behaves more like a pure force source with only damping, (that 
is, not necessarily resisting length changes to external force perturbations). Thus, 
even though there may still be force being generated during the movement (e.g., to 
oppose the external damping load in the experiments) the stiffness should still drop. 
Of course, these axe rough characterizations of the true complex behavior of muscle. 
Many of the finer details in the data are still hard to understand. For example, why 
do the stiffness changes slightly lead the muscle force changes (rising before the target 
is reached, that is, before peak acceleration)? This brings us to a consideration of the 
influence of reflexes. 

11.3.2 Reflex Contribution to Stiffness 

As the measured stiffnesses are low (giving natural frequencies below 3 Hz), it is 
possible that they are maintained by reflexes — in spite of feedback delays. That is, 
muscle forces during movement may be driven by reflexes to augment the effective 
mechanical stiffness of the arm (see Section 7.3). 

The contribution of reflexes can be assessed by measuring EMG activity (also see 
Section 7.6). An estimation of the transfer function between EMG and perturbation 
torques gives an indication of the reflex onset and magnitude. Fortunately, Soechting 
et al. [143] have measured such reflex transfer functions during ballistic forearm move¬ 
ments when applying PRBS force inputs with exactly the same frequency content as 
in this thesis (see Review Section 8.2). Their time-varying reflex impulse responses 
indicate that the reflex activity due to the perturbations is significant throughout the 
movement, with initial magnitudes and onset times similar to those of the same re¬ 
sponses measured while maintaining a fixed posture. The movement reflex responses 
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differ from the posture reflex responses in that they have significant activity at lags of 
up to 300 ms, indicating that proprioceptive information from as much as 300 ms in 
the past is used to control the movement — even information from before the onset 
of movement is used (re feedforward control). The reflex onset delays are of about 50 
ms, implying that the reflex loop gain can be significantly above unity for frequencies 
below 10 Hz (above 10 Hz the reflex gains must be below unity to ensure stability, 
see Section 1.1.1). Thus, it is possible that the reflexes can contribute to maintaining 
the stiffness values (<3 Hz natural frequency) measured in this thesis. 

11.4 Summary 

The contribution of this work has been to provide estimates of the time-varying me¬ 
chanical compliance of the arm during movement. The arm stiffness is modulated 
with the movement, dropping relatively low during motion. Such control of the me¬ 
chanical compliance could be advantageous from the standpoint of avoiding excessive 
contact forces in unexpected collisions, or more generally for producing constrained 
motions (see remarks in Chapter 12). 

The finding that the stiffness is low during movement is not inconsistent with an 
equilibrium trajectory hypothesis (i.e., static feedforward compensation). Simulations 
using the estimated stiffnesses show that, provided gravity is compensated for, the 
equilibrium trajectories reflect the actual trajectories closely (see Figure 7.15(b)). 
The reason for this is that, in spite of the equilibrium potential field being negligible 
during movement, the inertia of the arm continues the motion. 

The additional finding that dynamic scaling does not occur to compensate for a 
movement speed increase is inconsistent with the equilibrium trajectory hypothesis. 
McIntyre [162] anticipated this discrepancy, but argues that augmenting the equi¬ 
librium trajectory theory with a velocity referenced reflex feedback loop (including 
reasonable physiological delays) is sufficient to produce similar time scaled movements 
at different speeds. This solution preserves the simplicity of the theory, and prop¬ 
erly stresses the importance of velocity sensory information in the motor system (see 
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Chapter 12 


Future Research 

12.1 Kinematic Calibration 

12.1.1 Use of Force and Velocity Sensing 

In contrast to joint angle measurements, joint torque measurements enter linearly 
into the relation describing the endpoint force. Thus, kinematic calibration may be 
simplified by using force sensing. Specifically, it is possible to estimate the linearly 
occurring elements of the joint Jacobian matrix relating tip forces to joint torques. 
The kinematic parameters can then be analytically recovered from the Jacobian ma¬ 
trix elements. This analytic parameter estimation is possible because of the simple 
structure of the Jacobian matrix derived in Chapter 3. Similarly, as the joint an¬ 
gle Jacobian also relates joint velocities to endpoint velocities, it should likewise be 
possible to analytically estimate the kinematic parameters from velocity information. 

While analytic techniques may not give accurate parameter estimates when noise 
is present, they are important to develop, as a major difficulty in kinematic calibration 
is the estimation of initial parameter estimates. 

Finally, the use of velocity and force information in calibration is attractive from 
a biological point of view, as these quantities appear more readily available in the 
human motor system (re Chapter 6). 
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12.1.2 Global Uniqueness 

In Chapter 3 results pertaining to the uniqueness of the kinematic parameters were 
derived. It was shown that there are multiple solutions to the non-linear calibration 
equations, and that the equations become singular for certain parameter sets. As 
calibration must proceed iteratively, the following question should be addressed in the 
future. Is it possible to follow a continuous path in parameter space (by any iterative 
algorithm) from any initial point to a solution without crossing a singularity? It 
is reasonable to conjecture that the answer to this question is yes. At least for 
calibrating open kinematic chains, it is possible to show that the singularities divide 
up the parameter space into disjoint manifolds, each of which contains one of the 
2 n_1 solutions (provided the length parameters are suitably bounded). There may be 
additional solutions on each manifold, but there can only be a finite number of them. 
The manifolds provide a multiple covering of the output (endpoint) space. Thus, in 
the open-loop calibration case a solution can always be reached without crossing a 
singularity (i.e., without leaving a manifold). This needs to be elaborated and the 
full conjecture should be tested. 

12.1.3 Motor Psychophysics 

The relevance of model-based robotics control methods to human motor control re¬ 
quires further psychophysical experimentation. For example, experiments similar to 
the teleoperator adaption experiments described in Chapter 6 may help develop an 
understanding of how constrained our internal representations of sensorimotor trans¬ 
formations are. As well, simple experiments to characterize the accuracy of whole arm 
3-dimensional movement are necessary, before any strong conclusions can be drawn 
about how important calibration is to the motor system. 
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12.2 Airjet Experiments 


12.2.1 Compliance During Performance of Difficult Tasks 

Future airjet experiments should include characterizing the forearm’s compliance 
while performing more interesting tasks — such as catching a ball. The role of com¬ 
pliance in these tasks is more critical, and an understanding of how it is modulated 
may aid in building robots that can perform similar complex tasks. 

12.2.2 Role of Reflexes 

As mentioned in Chapter 7 the role of reflexes can be studied by estimating mod¬ 
els with sufficiently time lagged parameters — to account for delays. This can be 
done with the existing data. Alternatively, time-varying impulse responses can be 
estimated. As discussed in Chapter 9, this requires the development of an identifi¬ 
cation technique that incorporates a smoothness constraint on the impulse response. 
Finally, future experiments will record EMG during movement, and EMG-torque 
impulse responses will be identified. 

12.2.3 Whole Arm Compliance 

The main goal of the airjet project was (and still is) to characterize the 3-dimensional 
compliance of the arm during unconstrained movement. Generalization of the airjet 
perturbation device to a 3-dimensional version should be straight forward. Three- 
dimensional system identification techniques must also be developed. 

12.3 Role of Compliance in Constrained Arm Move¬ 
ments 

In previous work, not reported in this thesis, we performed a series of experiments 
to study the contact forces during constrained whole arm movement [163]. We found 
that a hypothesized strategy of relying on the compliance of the arm and ignoring 
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Appendix 


The damper torque load used in the airjet experiments has a torque-velocity char¬ 
acteristic given by Figure 12-1. In all cases the damper was set at 75 percent of 
maximum damping. 



Figure 12-1: Torque-velocity relation for the damper used to load the arm. The model 
used has the 500,000 cSt filling. At all times the viscosity was set at 75 percent of 
the maximum adjustable value. Damper manufactured by Kinetrol. 
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